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ABSTRACT 

A free  boundary  problem  is  a (steady-state)  boundary  value  problem 
involving  differential  equations  on  domains  parts  of  whose  boundaries,  the 
free  boundaries,  are  unknown  and  must  be  determined  as  part  of  the  solution. 
A trial  free-boundary  method  is  a numerical  method  for  solving  free  boundary 
problems  in  which  the  boundary  is  found  by  guessing  the  position  of  the 
boundary,  solving  a boundary- value  problem  in  the  resulting  region,  and 
using  this  solution  to  find  a new  guess  for  the  position  of  the  boundary,  the 


procedure  being  repeated  until  the  desired  accuracy  has  been  attained.  This 
report  gives  a comprehensive  survey  of  trial  free-boundary  methods. 

AMS  (MOS)  Subject  Classifications:  35J99;  35R99;  6502;  65N99 


Key  Words:  Free  boundary  problems,  Numerical  methods,  Trial  free- 
boundary  methods 
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A SURVEY  OE  TRIAL  FREE-BOUNDARY  METHODS  TOR  THE  NUMERIOAL 
SOLUTION  OF  FREE  BOUNDARY  PROBLEMS 
Colin  W.  Cryer 


0 . Introduction 

A free  boundary  problem  (EBP;  plural  FBPS)  is  a (steady-state) 
boundary  value  problem  involving  differential  equations  on  domains  parts 
of  whose  boundaries,  the  free  boundaries  (FB;  plural  FBS),  are  unknown  and 
must  be  determined  as  part  of  the  solution.  On  such  FBS  the  boundary 
conditions  needed  for  a fixed  boundary  value  problem  (a  problem  for  which 
the  boundaries  are  known)  are  supplemented  by  an  additional  boundary 
condition . 

The  problem  of  seepage  through  a rectangular  earth  dam  (see  Figure  1) 
illustrates  how  FBPS  arise. 

In  this  problem  water  from  an  upstream  reservoir  (or  head  water)  seeps 
through  a rectangular  earth  dam  to  a lower  reservoir  (or  tail  water).  The 
water  only  flows  through  the  region  £ and  the  remainder  of  the  dam 
remains  dry.  The  air/water  interface  is  denoted  by  r.  It  is  assumed 
that  the  dam  is  so  long  that  the  flow  is  two-dimensional. 

If  Darcy's  law  holds  then  the  flow  is  described  by  a velocity  potential 
6 which  is  related  to  the  hydraulic  head  h by 

* = -Kh  , (l) 
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Figure  1:  Seepage  through  a simple  rectangular  dam 


where  the  constant  K is  the  hydraulic  conductivity.  The  hydraulic  head 
satisfies  Laplace's  equation 


h + h =0,  in  & . 
xx  yy 

The  boundary  conditions  are: 

h = H,  on  AB  (interface  with  water  at  rest)  , 

= 0,  on  BC  (impervious  boundary)  , 

3n 

h = h ,,  on  CD  (interface  with  water  at  rest)  , 
d 

h = y,  on  DE  (interface  with  air)  , 

■^  = 0,  on  EA  (streamline)  . 

an  ’ 


(2) 


(3) 
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If  1'  wore  known  then  equation  (2)  with  boundary  conditions  (3) 
would  suffice  to  determine  h.  Since  r is  not  known  we  have  a FBP 
and  an  extra  condition  is  required  on  r.  This  extra  condition  is: 
h ~ y,  on  EA  (interface  with  air)  . 

As  illustrated  by  the  above  example,  a FBP  has  the  following  structure: 


£7u  = 0,  in  A , 

Bu  = 0,  on  ajp  , (5) 

Cu  = 0,  on  r f 

where  u is  the  solution  and  where  a represents  the  linear  or  nonlinear 

elliptic  differential  equation(s)  which  must  be  satisfied  by  u on  the 

open  set  t . B represents  the  linear  or  nonlinear  boundary  conditions  which 

mu^t  be  satisfied  on  the  boundary  9A  of  A.  If  9A  were  completely 

known  then  the  first  two  equations  of  (5)  would  define  u completely. 

However,  part  of  8A,  namely  the  FB  r,  is  unknown  and  C represents 

the  additional  linear  or  nonlinear  boundary  condition  imposed  on  r. 

In  (5)  it  is  to  be  understood  that  the  open  set  a may  consist  of 

m 

the  union  of  several  domains,  A = U a.,  in  which  case  the  first  equation 

i = l 1 

of  (5)  is  to  be  interpreted  as 


<7.u  = o in  A.,  1 < i < m . 
1 1 — — 


The  operators  £7.  may  or  may  not  be  the  same. 

The  domains  A.  will  be  called  free  domains  (FD;  plural,  FDS). 
When  it  is  necessary  the  number  oi  domains  will  be  indicated  by  saying 


that  the  problem  is  an  m- domain  FBP.  When  m = 1 we  sot  O E O , A E A 
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Using  this  terminology,  wc  may  describe  the  example  of  Figure  1 
as  a 1 -domain  FBI*  with  FD  and  FB  r. 


Most  of  the  FBPS  which  have  been  considered  in  the  literature  are 
1 -domain  FBPS.  A two-domain  FBP  would  arise,  for  example,  in  a coastal 
region  in  winch  porous  flow  of  fresh  water  occurred  in  a region  « and 
porous  flow  of  salt  water  occurred  in  a region  the  FB  bein9  the  salt 

water/fresh  water  interface. 

It  is  usually  neither  possible  nor  necessary  to  classify  the  boundary 
conditions  on  the  TB  into  fixed  boundary  conditions''  and  ' supple- 
mentary boundary  conditions",  and  the  division  of  the  boundary  conditions 
on  I'  into  those  represented  by  B and  that  represented  by  C is  an 
arbitrary  decision. 

The  FB  r is  defined  to  be  that  part  of  3£  whose  "shape"  is 
unknown;  the  remainder  of  3 £ will  be  called  the  fixed  boundary.  Points  at 
which  the  FB  intersects  the  fixed  boundary  will  be  called  points  of  detach- 
ment. If  the  location  of  a point  of  detachment  is  prescribed  it  is  said 
to  be  fixed;  otherwise  it  is  said  to  be  free.  We  illustrate  these  definitions 
by  considering  Figure  1 • The  point  A is  known  so  that 
A is  a fixed  point  of  detachment.  The  point  E is  not  known  so  that  E 

is  a free  point  of  detachment.  The  FB  is  EA;  although  _the  length  of DE  _js 

not  known  its  shape  is  known  so  that  DE  is  part  of  the  fixed  boundary 
and  not  part  of  the  FB . 


FBPS  have  recently  been  considered  in  control  theory,  and  this 
development  is  very  important  because  it  introduces  new  and  exciting 
ideas.  However,  the  vast  majority  of  FBPS  arise  in  continuum  mechanics. 

In  continuum  mechanics  a FBP  arises  in  principle  whenever  there  is  a 
line  or  surface  on  which  the  solution  is  discontinuous,  and  the  FB 
corresponds  to  this  line  of  discontinuity  or  surface  of  discontinuity  or 
singular  surface.  The  discontinuities  can  arise  for  many  reasons  including: 

(i)  Two  different  materials  may  be  present:  for  example, 
air  and  water. 

(ii)  The  material  may  occur  in  different  phases:  for  example, 
water  and  ice. 

(iii)  The  material  properties  may  depend  discontinuously  upon 
certain  variables : for  example,  the  material  may  change 
from  elastic  to  plastic  when  the  stress  reaches  a critical 
value. 

(iv)  There  may  be  a jump  in  the  pressure  or  velocity. 

FBPS  arise  in  the  following  branches  of  continuum  mechanics:  fluid 
mechanics,  porous  flow,  mechanics  of  solids,  heat  conduction  and 
diffusion,  electromagnetism,  and  gravitation. 

We  are  in  the  process  of  publishing  a series  of  surveys  on  the 
various  types  of  FBPS  (the  survey  on  porous  flow  FBPS  has  already 
appeared  (Cryer  [ 1976])).  We  are  also  in  the  process  of  publishing  a 
series  of  surveys  on  the  various  mathematical  and  numerical  aspects  of  FBPS. 
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The  present  survey  considers  the  class  of  numerical  methods  known 
trial  free  boundary  methods,  which  are  described  in  the  following 


section . 


1 . General  description  of  trial  tree  boundary  methods 
We  repeat  the  basic  FBP: 


d u - 0,  in  & 


/7u  = 0 , on  0$  , 


(1) 


C u = 0 , on  r . 

One  approach  to  solving  this  FBP  numerically  is  to  generate  a sequence  ot 

(k)  'k' 

trial  FBS  r and  approximations  u^/  to  the  corresponding  trial  solutions 

h 

( k) 

u as  follows: 


Step  0 . Guess  an  initial  trial  FB,  r 
Step  1 . 


(o) 


say. 


(k)  (k) 

Given  r let  & " be  the  corresponding  domain.  Compute 


(k)  (k) 

an  approximation,  u^  say,  to  the  solution  u of  the  problem 


(k) 


- (k> 

= o, 

Jk) 

du 

in  & , 

(2) 

&u(k) 

= 0, 

on  9*(k)  . 

* and 

(k) 

uh 

compute  a new  trial 

FBr(k+1)  by 

(k) 


requiring  that  C should  be  approximately  equal  to  zero  on 

r (k-i  1)  . „ „ , -(k)  . (k+1) 

r ; i.e.  move  the  boundary  from  r to  r 
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Following  Birkhoff  [ 19  61  ] we  will  call  this  method,  and  variations 
thereof,  a trial  tree  boundary  method-  (Previously,  we  have  called  methods 
of  this  type  "move-the-boundary"  methods  (Cryer  [ 1968,  1970a])  but 
Birkhoff's  terminology  is  more  descriptive.) 

, 00 

In  Step  Z the  computation  of  the  approximate  solution  u requires 

the  approximate  solution  of  a fixed  boundary  value  problem,  which  can  be 

done  using  any  standard  method  such  as  finite  elements  or  finite 

differences.  We  then  speak  of  a trial  free  boundary  method  using  finite 

elements,  a trial  free  boundary  method  using  finite  differences,  etc. 

Steps  1 and  2 in  a trial  free  boundary  method,  namely  the  computation 

(k) 

of  the  approximate  solution  u^  and  the  generation  of  the  new  trial 
FB  are  essentially  independent  and  are  considered  separately 

in  sections  Z and  3,  respectively.  In  section  3.4  we  summarize  the 
numerical  experiences  of  workers  who  have  used  trial  free  boundary 
methods  and  make  some  suggestions  regarding  their  use.  Finally,  error 
estimates  and  related  questions  are  considered  in  section  4. 

The  advantages  of  trial  free  boundary  methods  are: 

(i)  Trial  free  boundary  methods  are,  in  principle,  applicable  to  all 
FBPS  and  require  no  preliminary  analysis. 

(ii)  The  main  computational  effort  is  in  Step  1 which  involves  the 
solution  of  a linear  problem  if  Q is  linear. 

(iii)  All  computations  are  carried  out  in  the  physical  domain  so  that  "physical 
intuition"  can  be  used  and  a "feeling"  for  the  solution  can  develop 


as  the  computation  proceeds. 
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Because  or  these  advantages,  trial  tree  boundary  methods  are  particularly 
well  suited  for  hand  computation,  an  1 before  the  advent  of  computers 
trial  tree  boundary  methods  were  the  most  widely  used  numerical  methods 
for  solving  FBPS. 

Unfortunately,  trial  free  boundary  methods  also  have  several 
disadvantages: 


(a) 

(b) 


To  implement  Step  1 one  must  be  able  to  compute  the  approximate 

solutions  uUj  for  the  general  domains 
h 

It  is  difficult  to  formulate  Step  2 precisely.  Each  problem  seems 


to  require  a different  technique  and  while  an  individual  research 


worker  can  accumulate  experience,  it  is  hard  to  transmit  this 
knowledge  to  others.  Many  cases  of  non-convergence  have  been 
reported. 

(c)  It  is  hard  to  achieve  high  accuracy  and  it  is  hard  to  estimate 
(k)  (k) 

the  error  in  u^  and  r . In  some  problems  the  shape  of  the 
FB  is  very  sensitive  to  small  errors  in  the  condition  Cu  = 0,  and 
it  is  hard  to  achieve  high  accuracy  near  points  of  separation.  In 
general,  it  is  difficult  to  obtain  accuracy  greater  than  about  1%  . 
Because  of  these  disadvantages  it  proved  relatively  difficult  to  implement 
trial  free  boundary  methods  on  computers,  particularly  on  computers  with 
small  memories.  There  were  a few  successful  implementations  during 
the  early  1960's,  but  it  is  only  since  about  1967  that  there  has  been  a 
substantial  number  of  successful  implementations. 


for  a number  of  years  we  felt  that  trial  free  boundary  methods 


would  and  should  be  replaced  by  alternative  methods  such  as  numerical 
methods  based  upon  variational  inequalities.  However,  we  have  recently 
changed  our  opinion  for  the  following  reasons: 

(а)  There  are  many  FBPS,  such  as  viscous  flow  fluid  mechanics 
FBPS,  for  which  trial  free  boundary  methods  are  at  present  the 
only  available  method. 

(3)  The  difficulty  of  implementing  Step  1 has  been  reduced  with 
the  use  of  finite  elements  and  other  techniques  which  are 
more  flexible  than  finite  differences. 

(y)  Step  2 has  recently  been  successfully  implemented  in  several 
applications  using  "integral  methods"  (see  section  3.2.2). 

(б)  Although  error  estimates  and  convergence  proofs  are  still 
lacking,  there  have  been  a large  number  of  successful 
applications . 

We,  therefore,  believe  that  trial  free  boundary  methods  will  remain  a 
useful  and  competitive  approach  to  the  numerical  solution  of  certain 
classes  of  FBPS. 
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. A 


: . Step  1 : Computation  of  the  approximate  trial  solution 

( k ) 

The  approximate  trial  solution  u,  ' is  the  approximate  solution  of 

h 


the  fixed  boundary  value  problem 

,(k) 


t7u 


= 0,  in  fi 


(k) 


Bu.  = 0,  on 


(k) 


(1) 


.(k) 


and  u.  ' can  be  computed  by  any  approximation  method  for  fixed  boundary 
n 

value  problems. 

The  term  approximation  method  includes  both  numerical  methods 
(such  as  the  method  of  finite  differences)  and  nonnumerical  methods  (such 
as  analog  methods).  Each  of  the  following  subsections  discusses  the 

, k) 

solution  of  u,  by  a particular  approximation  method, 
h 

Birkhoff  [ 1971  ] gives  a very  readable  introduction  to  numerical 
methods  for  fixed  boundary  value  problems.  Extensive  bibliographies  (not 
restricted  to  elliptic  equations)  are  given  by  Voight  [ 1967  ] , Giese  [ 1969  ] , 
and  Forsythe  [ 1971  J;  the  bibliography  of  Giese  provides  an  exhaustive 
coverage  of  early  work.  Among  the  many  symposia  on  the  subject  are 
Birkhoff  and  Varga  [ 1969  ),  and  Hubbard  [ 1976]  . Specialized  texts  discussing 
particular  approximation  methods  are  referred  to  in  the  appropriate 


subsections  below. 

In  the  remainder  of  this  section  we  make  some  general  remarks  about 

(k) 

the  relative  efficiency  of  different  methods  of  computing  u^  . The  relative 
efficiency  of  approximation  methods  for  FBPS  is,  of  course,  closely  related  to 
their  relative  efficiency  for  fixed  boundary  problems,  but  there  are  some 
interesting  differences  which  deserve  further  study: 
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(i)  There  is  growing  interest  in  comparing  the  efficiency  of  numerical 


w 


methods  for  solving  the  fixed  boundary  value  problem  (l). 

Terry  [ 1967  ] has  compared  integral  equation  methods  and  finite 

difference  methods  for  Neumann  boundary  value  problems  and  finds 

that  the  finite  difference  methods  are  faster.  It  is  very  likely 

that  comparisons  between  finite  difference  methods  and  finite  element 

methods  for  solving  fixed  boundary  value  problems  will  soon  be  available. 

(ii)  Ail  the  numerical  methods  for  solving  (1)  require  the  solution  of  a 

system  of  n algebraic  equations  of  the  form 

(h)  = B(k)  (2) 

h h 

i k)  ( k)  n 

where  U and  B are  in  n-dimensional  Euclidean  space  R and 
h h 

A) : Rn  - Rn.  The  equations  (2)  are  linear  if  O and  B in  (1)  are 
n 

linear;  otherwise  the  equations  (2)  are  nonlinear. 


Basic  references  on  the  solution  of  systems  of  linear  algebraic 

equations  include:  Varga  [ 1962],  Wilkinson  [ 1965],  Young  [ 1971  ],  and 

Stewart  [ 197  3].  Basic  references  on  the  solution  of  systems  of  nonlinear 

algebraic  equations  include:  Ortega  and  Rheinholdt  [ 1970],  Rabinowitz  [ 1970], 

Byrne  and  Hall  [ 197  3],  and  Rheinboldt  [ 1974  ] . 

(k) 

If  A^  is  linear  then  the  numerical  methods  for  solving  (2)  fall 
into  two  categories:  direct  methods  and  iterative  methods.  The  relative 
efficiency  of  direct  methods  and  iterative  methods  for  solving  (2)  depends 
upon  the  structure  of  A]  and  computing  facilities  which  are  available. 


improved.  For  example,  if  the  equations  are  suitably  ordered,  then  many 

(k) 

of  the  computations  required  to  solve  for  U using  Gaussian  elimination 

h 

( J^—  1 ) 

will  already  have  been  performed  when  solving  for  U'  , so  that 

h 

judicious  preservation  of  previous  results  should  lead  to  a substantial 

reduction  in  computing  time. 

(iii)  There  is  a close  relationship  between  the  popularity  (and  efficiency  ?) 
of  a given  approximation  method  for  solving  (l)  and  the  computing 
facilities  available.  For  example,  integral  equation  methods  and 
finite  element  methods  only  became  popular  when  computers  with 
large  high  speed  memories  became  generally  available.  This  rela- 
tionship between  approximation  methods  and  computing  facilities 
is  clearly  apparent  from  the  references  quoted  in  this  survey 
which  span  a period  of  sixty  years. 
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The  relationship  between  numerical  methods  and  computing 


facilities  has  not  been  very  obvious  in  recent  years  because  during  the 
period  1960-1975  the  computers  which  have  been  available  have  been 
rather  similar  from  a numerical  (as  opposed  to  systems)  viewpoint.  A 
typical  computer  has  had  a microsecond  arithmetic  unit,  a reasonably 
large  microsecond  core  store,  and  a large  secondary  store.  A new  era 
is  now  beginning  in  which  a much  greater  variety  of  computers  will  be 
available.  Already,  the  parallel  computer  ILLIAC  IV  which  has  64 
separate  arithmetic  units  has  been  used  to  solve  FBPS  (MacCormack  and 
Stevens  [ 1974]).  A corresponding  re-evaluation  of  numerical  methods 
will  occur.  To  mention  but  one  possibility,  when  using  a trial  free 
boundary  method  it  may  be  most  efficient  to  make  several  guesses  for 
the  FB  and  solve  the  corresponding  fixed  boundary  value  problems 
simultaneously. 

It  is  impossible  to  predict  long-term  future  developments  in  computers, 

and  it  is  quite  possible  that  approximation  methods  which  are  at  present 

regarded  as  out-of-date  will  again  become  popular.  To  mention  one 

possibility,  suppose  that  it  became  possible  to  build  a picosecond 

(10  second)  computer  but  that  because  of  technical  and  financial 

restrictions  the  storage  consisted  of  a l 000-word  picosecond  store  with 

-9 

a million-word  nanosecond  (10  second)  secondary  store.  This  computer 
would  be  109  times  faster  than  the  computers  of  19  50,  but  the  storage  alloca- 
tion problems  would  be  very  similar.  Iterative  methods  of  solving  (2),  which 
require  less  storage  than  direct  methods,  would  again  become  very  popular. 
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(iv)  Most  approximation  methods  for  fixed  boundary  value  problems 


are  based  on  the  assumption  that  the  solution  is  smooth.  This 
assumption  is  usually  true  if  the  boundary  is  smooth,  but  most 
FBPS  involve  boundaries  with  corners  in  the  neighbourhood  of 
which  the  solution  can  have  singularities.  If  the  corner  is  a point 
of  intersection  of  the  FB  with  the  fixed  boundary,  then  the  analysis 
of  the  singularity  is  even  more  difficult. 

In  the  literature  on  trial  free  boundary  methods  there  has  been  little 
discussion  of  the  treatment  of  corner  singularities.  This  is  clearly  a 
subject  which  requires  further  study,  and  we  discuss  it  briefly  in  section  4.2. 
(v)  All  numerical  methods  for  the  fixed  boundary  value  problem  (l) 
require  that  a decision  be  made  about  the  location  of  certain 
special  points  which  we  will  call  "meshpoints ".  Examples  of 
meshpoints  are:  netpoints  (for  finite  difference  methods);  nodes 
(for  finite  element  methods);  quadrature  points  (for  integral  equation 
methods);  and  collocation  points. 

It  is  desirable  that  the  choice  of  meshpoints  should  be  made  auto- 
matically by  the  computer  rather  than  by  hand  for  several  reasons: 

(a)  When  using  a trial-free-boundary  method,  the  domain  changes 
at  every  iteration. 

(b)  If  the  solution  is  ill-behaved  in  a certain  subregion  of  fl, 
as  is  often  that  case  for  FBPS,  it  is  desirable  to  be  able 
to  automatically  refine  the  mesh  in  this  subregion. 
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(c)  The  most  effective  check  of  numerical  results  is  to  repeat 
the  computation  with  a finer  mesh.  The  labor  involved  in 
data  preparation  increases  rapidly  as  the  mesh  is  refined, 
and  if  the  data  is  prepared  by  hand  the  necessary  checks 
may  not  be  carried  out. 

Although  automatic  mesh-generation  programs  are  feasible,  almost  all 

programs  for  trial  free  boundary  methods  have  used  a substantial  amount 

of  manually  prepared  data.  It  is  to  be  hoped  that  this  will  soon  be  remedied. 

It  is  also  desirable  that  the  meshpoints  be  chosen  so  that  the 

approximations  u^  change  "continuously"  as  k increases.  For 

example,  Arms  and  Gates  [ 1957]  used  a trial  free  boundary  method  with 

finite  differences  and  the  same  grid  for  all  k.  Arms  and  Gates  [ 19  57,  p.  6] 

.(k-1)  (k) 

observe  that  discontinuities  occurred  whenever  "moving"  r to  r 

involved  crossing  a gridpoint,  and  these  discontinuities  invitiated  their 
method . 

(vi)  Finally,  when  uj^  has  been  computed  it  is  necessary  in  Step  2 

to  compute  Cuj^ , and  it  is  of  course  desirable  that  this  computa- 
tion should  be  relatively  easy. 

As  is  clear  from  the  above  remarks,  it  is  not  possible  to  indicate 

which  approximation  method  is  "best"  for  Step  1 . For  several  years 

finite  differences  were  widely  used.  However,  it  is  relatively  difficult 

( J^) 

to  implement  finite  difference  methods  to  handle  a general  domain  8 , 

and  it  is  not  surprising  that  finite  element  methods,  which  are  much  more 
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flexible,  have  recently  become  very  popular.  Integral  equation  methods 


and  i allocation  methods  are  also  very  flexible,  and  we  anticipate  that 
their  use  will  increase. 

2 . 1 The  trial  free  boundary  method  with  finite  differences 

finite  difference  methods  for  elliptic  equations  were  introduced 
by  Runge  in  1908  and  Richardson  in  1910.  One  would  have  expected  that 
finite  difference  methods  would  have  been  quickly  applied  to  FBPS, 
particularly  since  Runge  was  familiar  with  the  idea  of  trial  free  boundary 
methods  for  FBPS.  (Blasius  [ 1 9 1 0 , p.  97]  acknowledges  a suggestion  made 
by  Runge  in  the  Gottingen  seminar  on  hydrodynamics.)  We  have  consulted 
many  of  the  early  papers  of  finite  difference  methods  (such  as  those 
listed  by  Giese  [ 1 9 69  ] ) , but  the  first  application  of  trial  free  boundary 
methods  using  finite  differences  is  apparently  due  to  Vandrey  [ 1940]  who 
computed  the  flow  through  a circular  Borda  mouthpiece. 

Independently,  Shaw  and  Southwell  [ 1941  ] used  finite  differences 
to  solve  the  problem  of  seepage  through  an  earth  dam  for  three  geometries. 
Southwell  and  Vaisey  [ 1946]  applied  this  method  to  eleven  fluid  mechanics  FBPS. 
The  work  of  Southwell  and  his  colleagues  was  followed  by  many  similar 
computations  by  other  workers.  Not  every  author  indicates  how  much 
effort  was  required,  but  some  do:  Yang  [ 1949,  p.  54],  who  solved  two 
porous  flow  FBPS,  states  that  his  first  problem  took  224  hours  and  that 
subsequent  problems  took  about  60  hours;  Rouse  and  Abul-Fetouh  [ 19  50,  p.  42  3] 


W 


who  solved  the  axisymmetric  wall  jet  FBP,  state  that  their  first  problem 
took  l SO  hours.  We  should  be  duly  grateful  for  computers. 

The  implementation  of  trial  free  boundary  methods  on  early  computers 
was  not  a straightforward  matter  for  two  reasons:  (i)  early  computers 
had  smail  memories  which  made  the  implementation  of  Step  1 difficult; 

(ii)  tt  was  necessary  to  automate  Step  2,  and  it  was  far  from  clear  how 
this  should  be  done. 

The  first  attempt  to  use  a digital  computer  is  apparently  due  to 
Young  et  al.  [ 1 9 55  ] and  Arms  and  Gates  [ 1 9 57  ] . This  work  was  not 
entirely  successful.  Jakobsson  and  Floberg  [ 1957  ] solved  several  lubrica- 
tion cavitation  problems;  they  used  a computer  but  moved  the  FB  by  hand 
(see  Jakobsson  and  Floberg  [ 1957,  p.  33]).  In  the  applications  listed  at 
the  end  of  this  section,  all  the  results  obtained  after  1 9 6 1 were  obtained 
using  a computer. 

It  is  appropriate  to  mention  here  that  we  were  once  told  that  an 
IMB  650  program,  Program  4.0.010  of  February  1962,  was  capable  of 
solving  porous  flow  FBPS  using  trial  free  boundary  methods  with  finite 
differences,  but  we  have  been  unable  to  locate  this  program. 

The  basic  texts  on  finite  difference  methods  are  the  books  of  Forsythe 
and  Wasow  [ 1960]  and  Coilatz  [ I960  ] . Greenspan  [ 1968]  discusses  several 
applications . Remson,  Hornberger,  and  Molz  [ 1971  ] give  applications  to 
porous  flow. 
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The  finite  difference  method  consists,  in  essence,  of  the  following 


steps: 

no  (k) 

(i)  Choose  a set  V of  n netpoints  in  O U 9fl'  . 

(k) 

(ii)  Set  up  a system  of  n equations  connecting  the  values  U, 

h 

, (k) 

of  the  approximate  solution  u^  at  the  netpoints, 

A<kVk)  = B<k)  . 

h h h 

(iii)  Solve  this  system  of  n equations. 

The  netpoints  are  usually  taken  to  be  the  intersections  of  an  array 

(k) 

of  gridlines  with  each  other  and  with  9jfl  . The  implementation  is 

simplified  if  the  network  is  gridlike:  that  is,  the  gridlines  intersect  the 

boundary  9£  only  at  gridpoints  (intersections  of  gridlines).  The  advantage 

of  a gridlike  network  is  that  there  are  no  awkward  "special  netpoints" 

near  the  boundary  which  require  special  treatment. 

If  the  gridlines  are  equally. spaced  then  it  is  obviously  often  not 

(k)  P (k) 

possible  to  choose  the  network  so  that  it  is  gridlike.  If  9fr  = U 9.S 

J = 1 3 

(k)  (k)  (k) 

where  dfi  is  a curve  and  9„fi  B & are  line  segments,  one 

1 2 p 

(k) 

can  Jioose  a number  of  points  on  9^£  and  define  the  gridlines  to 
be  the  (unequally  spaced)  horizontal  and  vertical  lines  passing  through 
these  points;  with  luck,  the  resulting  network  will  be  gridlike.  Figure  1 
shows  such  a grid  for  the  problem  of  a jet  issuing  from  a container 
(Cryer  [ 1968]).  However,  Cryer  [ 1969,  1971]  has  given  examples  of 
quite  simple  domains  ^ for  which  no  gridlike  network  exists. 


1) 
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Figure  1 : Gridlines  defined  by  boundary  points 


Handworkers  commonly  used  a graded  grid,  that  is,  a grid  in  which 
the  basic  grid  consisted  of  equally  spaced  gridlines  with  gridlength  h, 
but  in  which  smaller  gridlengths  were  used  near  comers  and  other 
singular  points.  Figure  2 (Southwell  and  Vaisey  [ 1946,  p.  140],  Southwell 
[ 1946,  Figure  11  6])  shows  part  of  a graded  grid  for  flow  from  an  orifice. 

By  hand  it  is  easier  to  use  a graded  grid  (as  in  Figure  2)  rather  than 
an  unequally  spaced  grid  (as  in  Figure  l)  because  the  finite  difference 
formulas  are  simpler.  On  a computer,  it  is  easier  to  use  an  unequally 
spaced  grid  rather  than  a graded  grid  because  the  computer  can  easily 
handle  the  more  complicated  finite  difference  formulas  for  an  unequally 
spaced  grid  but  encounters  storage  allocation  problems  with  a graded  grid. 


-20- 


So  far  as  we  are  aware,  only  Rippln  [ 19  99]  and  Rippin  and  Davidson  | 1967) 


have  used  a graded  grid  in  a computer  implementation  of  a trial  free 
boundary  method. 

It  would  of  course  be  possible  to  use  a general  triangular  mesh 

i Kellogg  [ 19  6-4))  but  then  one  might  as  well  use  a finite  element  method. 

The  development  of  programs  which  can  automatically  generate  a 

mesh  is  laborious  but  feasible.  The  author  (Cryer  [ 1962,  1968,  l 9 7 0 a J ) 

developed  such  a program  for  use  with  trial  free  boundary  methods. 

00 

The  matrix  A,  in  (1)  is  sparse  and  has  a great  deal  of  structure, 
h 

There  are  many  possible  methods  of  solution  which  take  advantage  of 
this  structure.  The  main  methods  are: 

( a)  Iterative  methods 

(a)  Classical  pointwise  methods  such  as  SOR  (systematic  over 
relaxation).  See  Varga  [1962],  Young  [1971]. 

(p)  Splitting  methods.  See  Yanenko  [1971],  Stone  [1968]. 

(b)  Direct  methods 

(a)  Fast  direct  methods . See  Dorr  [ 1970  ] , Buzbee  and  Dorr  [ 1974  ] . 

((3)  Sparse  matrix  methods.  See  Willoughby  [ 1968],  Reid  [1971], 

Tewerson  [ 197  3 ] , Bunch  and  Rose  [ 197  6 ] . 

As  discussed  in  section  2,  for  each  method,  whether  iterative  or 

direct,  the  question  arises  as  to  whether  the  method  can  be  adapted  to  take 

1 

advantage  of  the  iterative  nature  of  trial  free  boundary  methods. 
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The*  following  TBPS  have  been  solved  using  the  trial  free  boundary 


method  in  conjunction  with  finite  differences: 

I'll  lid  me  chanics 

Plane  and  axially  symmetric  cusped  cavities:  Southwell  and  Vaisey  f 1946]. 

Plane  Riabouchinsky  cavity:  Mogal  and  Street  | 197  2,  19’.  4]  . 

Axially  symmetric  Riabouchinsky  cavity:  Brunauer  [1951]  , Young,  Gates, 
Arms,  and  Eliezer  [ 19  55]  , Arms  and  Gates  [ 19  57]  , Sankar  [ 1967]  , 
Fox  and  Sankar [ 197  3] . 

Plane  Borda  mouthpiece:  Southwell  and  Vaisey  [ 1946]  . 

Axially  symmetric  Borda  mouthpiece:  Vandrey  [ 1940]  , Southwell  and 
Vaisey  [ 1946]  . 

Plane  jet  from  an  orifice  in  a tube:  Cryer  [ 1968]  . 

Axially  symmetric  jet  from  an  orifice  in  a tube  (with  and  without  gravity): 
Southwell  and  Vaisey  [ 1946],  Abul-Fetouh  [ 1949],  Rouse  and 
Abul-Fetouh  [ 19  50]  . 

Flow  over  a plane  weir:  Southwell  and  Vaisey  [1946]  . 

Flow  over  a circular  weir:  Citrini  [19  51]  . 

Flow  over  a submerged  plate:  Dumitrescu,  Ionescu  and  Toth  [ 1956]  . 

Flow  over  a submerged  arc:  Dumitrescu,  Ionescu,  and  Toth  [ 1956]  . 

Flow  over  a submerged  vortex:  von  Kerczek  [ 196  5]  . 

Flow  under  a planing  surface:  Southwell  and  Vaisey  [ 1946]  , von  Kerczek 
[ 1965]  (unsuccessful). 


How  under  a sluice  gate:  Southwell  and  Vaiscy  [ 1946]  . 

Plano  confined  infinite  symmetric  bubbles:  Stewart  ( 1965]  , Stewart  arid 


Davidson  f 1967  ] . 

Axially  symmetric  confined  infinite  bubbles : Duinitrescu,  Ionescu,  and 
Toth  [ 19  56]  , Stewart  [ 196  5]  , Stewart  and  Davidson  [ 1967]  . 

Axially  symmetric  confined  spherical  cap  bubbles  with  wake:  Rippin  [ 1959] 
and  Rippin  and  Davidson  [ 1967]  . 

Periodic  progressing  gravity  waves:  Southwell  and  Vaisey  [ 1946]  . 

Solitary  waves:  Southwell  and  Vaisey  [ 1946]  . 

Two  and  three-dimensional  aerofoils:  Allen  and  Dennis  [ 19  53]  . 

Swirling  flow  through  an  orifice:  Binnie  and  Davidson  [ 1949]  . 

Flame  in  a channel:  Ball  [ 19  51]  . 

Lubrication  cavitation:  Jakobsson  and  Floberg  [ 19  57 ] . 

Porous  flow 

Simple  rectangular  dam:  Shaw  and  Southwell  [1941,  p.  9],  McNown, 

Hsu,  and  Yih  [ 19  55,  p.  668],  Finnemore  and  Perry  [ 1968  , p.  1066], 
J.  M.  Taylor  [1971,  p.  56]  and  Aitchison  [ 197  2]  , Comincioli, 
Guerri,  and  Volpi  [ 1971]  . 

Simple  trapezoidal  dam:  Yang  [ 1949,  p.  48]. 

Polygonal  dams:  Shaw  and  Southwell  [1941]  and  Murray  [i960,  p.  156] 
(trapezoidal  dams  with  toe  drains  on  pervious  sublayers); 

Freeze  [ 1971a]  (inhomogeneous  trapezoidal  dams  for  saturated- 
unsaturated  flow). 

Drainagei  van  Dcemter  [ 19 50]  . 
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Magnetohydro  Jynamics 

Stellar  evolution:  Cryer  [ 1962,  1 9 68  ] . 

Tokamak  and  toroidal  compressor:  Stevens  [ 1974]. 

Electrohydrostatics 

Fluid  film  in  an  electric  field:  Michael  and  O'Neill  [ 1972a]. 

2.2.  The  trial  free  boundary  method  with  finite  elements 

The  finite  element  method  is  usually  attributed  to  Courant  [ 1943] 
although  Le  Roux  [ 1914]  derived  the  finite  difference  equations  for 
Laplace's  equation  by  minimizing  the  Dirichlet  integral  over  piecewise 
linear  functions.  The  first  serious  computations  were  performed  by 
Turner,  Clough,  Martin,  and  Topp  [ 19  56]  and  Clough  [i960].  The  first 
application  of  finite  elements  to  trial  free  boundary  methods  is  probably 
due  to  R.  L.  Taylor  [ 1966]  who  solved  several  porous  flow  FBPS. 

Information  about  the  method  of  finite  elements  is  given  by  Birkhoff 
[ 1971  ] , Varga  [ 1 97 1 ] , Zienkiewicz  [ 1971  ] , Aubin  [ 1972],  Aziz  [ 1972  | , 
Desai  and  Abel  [ 1972],  Oden  [ 1972],  Whiteman  [ 1 97  3,  1975],  Strang  and 
Fix  [ 1973],  de  Boor  [ 1974],  Ciarlet  [ to  appear] . 
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Because  of  their  greater  flexibility,  finite  element  methods  encounter 


fewer  topologi  al  difficulties  than  finite  difference  methods  when  used 
with  trial  free  boundary  methods.  It  is  of  course  desirable  to  choose 
the  finite  elements  so  as  to  take  account  of  any  singularities.  In  addition 
to  the  discussions  in  textbooks  ( Zienkiewicz  [ 1 97  l ) , Desai  and  Abel  [ 1972  ]), 
the  paper  of  Iversen  [ 1969)  gives  a clear  discussion  of  the  problems 
involved. 

When  using  finite  elements  in  conjunction  with  a trial  free  boundary 
method  it  is  customary  to  choose  the  finite  element  array  a-priori.  The 
same  finite  array  is  used  throughout  the  computation:  as  the  trial  free 
boundary  changes  only  the  coordinates  of  the  finite  elements  adjacent  to 
the  FB  change.  Of  course,  if  the  trial  free  boundary  changes  a great  deal, 
then  an  element  may  become  very  distorted  or  may  even  disappear,  and 
this  must  be  checked  for.  A typical  finite  element  array  is  shown  in 
Figure  1 for  the  problem  of  viscous  flow  from  an  axially  symmetric  pipe 
(Nickell,  Tanner  and  Caswell  [ 1974,  Figure  5a]);  the  finite  elements  have 
been  chosen  so  as  to  take  account  of  the  singularity  at  the  point  A where 
the  fluid  leaves  the  pipe. 

There  are  70  elements  in  Figure  1 . Nickell,  Tanner  and  Caswell  [ 1974] 
do  not  indicate  how  many  equations  were  involved,  but  in  a related  problem 
(the  stick-slip  problem)  150  elements  involved  1000  equations,  so  that  we 
may  conclude  that  the  array  shown  in  Figure  1 involved  about  500  equations. 
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figure  1:  An  array  of  quadrilateral  elements 


It  is  typical  of  finite  element  solutions  of  FBPS  that  relatively  coarse 
arrays  are  used.  For  comparison,  the  finite  difference  grid  part  of  which 
is  shown  in  Figure  2 of  section  2.  1 involved  1028  equations  which  were 
solved  by  hand. 

Winslow  [ 19  67]  describes  a method  for  the  automatic  triangulation 
of  a region,  and  this  method  has  been  extended  by  Godunov  [ 1971]; 

George  [ 1971,  p.  27]  gives  a critical  discussion  of  Winslow's  method. 

George  [ 1971,  chapter  2]  has  an  excellent  discussion  of  the 
problems  arising  when  automatically  generating  a triangulation  of  a domain. 
George  reviews  previous  work  in  the  literature,  suggests  two  algorithms, 
and  provides  a program  listing. 

Sewell  [ 1972]  describes  an  algorithm  for  automatically  refining  a 
triangulation  of  a polygonal  region  in  R . The  algorithm  is  designed  so 
that:  (i)  no  triangle  A has  a very  small  angle;  (ii)  the  integral 
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- 2/(n+2) 

dx  dy  * 1 ) 

is  approximately  equal  for  each  triangle  A,  where  the  solution  u is 
being  approximated  by  polynomials  of  degree  n.  A sequence  of  triangula- 
tions *rr  is  generated:  the  triangulation  *k  is  used  to  compute  an 
approximation  u(k)  from  which  the  integrand  in  (1)  can  be  estimated 

prior  to  generation  of  the  refinement 

The  finite  element  method  leads  to  a system  of  linear  equations 


(k)  (k)  = (k)  (2) 

h h h 

r 

The  matrix  a!^  is  sparse  and  has  a band  structure,  but  has  in  general 
h 

much  less  structure  than  the  corresponding  matrix  for  the  finite  difference 
method. 

In  many  programs  the  system  (2)  is  solved  using  a band  matrix  version 
of  Gaussian  elimination.  A number  of  different  ideas  have  been  used: 

(a)  Instead  of  requiring  a band  of  uniform  width,  one  can  allow  bands 

of  variable  width.  This  is  a straightforward  generalization:  among  the 
first  implementations  was  that  of  Cryer  [ 1962,  Chapter  5]  . 

(b)  Instead  of  generating  all  the  equations  at  once,  one  can  use  the 
fmnt.l  elution  method  (Irons  [ 1970])  in  which  the  equations  are 


generated  when  needed  in  the  Gaussian  elimination. 

(c)  One  can  try  to  re-order  the  equations  so  as  to  reduce  the  band  width 
George  [197  2,  197  4]  has  suggested  some  interesting  ideas  here. 
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The  solution  by  direct  methods  of  sparse  linear  systems  is  at 

present  a very  active  field  of  research.  See  the  books  and  conferences: 

Willoughby  [ 19 68  | , Reid  [1971],  Tewerson  [ 197  3 ] , Bunch  and  Rose  [ 1976  J . 

In  a few  cases  the  system  (2)  has  been  solved  by  an  iterative 

methods.  For  trial  free  boundary  methods  this  is  a desirable  approach 

( k-1) 

because  the  previous  solution  is  available.  Unfortunately, 

( k) 

the  lack  of  detailed  structure  in  the  matrix  A'  makes  it  difficult  to 

h 

analyse.  For  example,  Fix  and  Larsen  [ 1971]  have  considered  S.O.R. 

methods  but  in  his  review  of  the  paper  Varga  (MR  4_5  (197  3)  ?29  35)  observes 

that  the  main  theorem  of  Fix  and  Larsen  contains  a nontrivial  flaw,  and 

adds  that  he  doubts  whether  the  convergence  properties  claimed  are  valid 

without  additional  hypotheses.  Schultchen  [ 197  3]  has  tested  an  immense 

( k. ) 

number  of  interative  methods  for  finite  element  matrices  A,  and 

h 

concludes  that  a particular  version  of  the  conjugate  gradient  method  is  the 
most  efficient  procedure.  See  also  Fried  [1970], 

It  is  also  important  to  have  information  about  the  round-off  errors 
which  arise  when  solving  the  system  (2).  See  Roy  [ 1971  ] , and  S.  K.  Yang  [1974 
The  following  FBPS  have  been  solved  using  the  trial  free  boundary 
method  in  conjunction  with  the  finite  element  method: 


Fluid  mechanics  FBPS 

General  references  on  the  application  of  finite  elements  to  fluid 
mechanics  problems  include:  de  Vries  and  Norrie  [ 197  3 ] , Oden,  Zienkiewicz, 
Gallagher,  and  Taylor  (1974). 
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The  following  specific  problems  have  been  solved: 

Cavitating  hydrofoils:  Dobray  and  Baker  [1974].  (This  reference  was 
kindly  drawn  to  our  attention  by  B.  E.  Larock.  ) 

Viscous  flow  from  a tube  (die-swell):  Tanner  [ 1973]  and  Nickell,  Tanner, 
and  Caswell  [1974].  (These  references  were  kindly  drawn  to  our 
attention  by  J.  Liu.) 

Plane  flow  from  symmetric  polygonal  nozzles:  S.  T.  K.  Chan  [1971]  and 
S.  T.  K.  Chan,  Larock  and  Herrman  [ 197  3]. 

Axisymmetric  flows  from  orifices  and  valves:  Chang  and  Larock  [ 1973]. 
Three-dimensional  jet  in  a transverse  gravity  field:  Larock  [to  appear]. 
Circular  crested  weir  (spillway):  Tu  [1971]. 

Gated  spillway  crest:  Larock  [to  appear]. 

Plane  jet  impinging  on  a plane:  Tu  [1971]. 

Surface  wave  generated  by  a vortex  (with  linearized  equations):  Tu  [1971]. 
Porous  flow  FBPS 

R.  L.  Taylor  [19  66]  gives  a general  program  for  solving  plane  and 
axisymmetric  porous  flow  FBPS  using  finite  elements.  A later  version  of 
this  program  is  given  by  Kealy  and  Busch  [1971]. 

Variational  principles  for  porous  flow  FBPS  have  been  developed  by 
Mauersberger  [ 1965,  1965b].  Ponter  [ 1972]  has  derived  dual  minimum 
principles  for  porous  flow  fixed  boundary  problems  and  used  them  in  finite 
element  computations;  presumably  these  dual  minimum  principles  could 
also  be  used  for  FBPS. 
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The  following  specific  problems  have  been  solved: 

Simple  rectangular  dam:  Kealy  and  Busch  [ 1971,  p.  9],  Kealy  and 
Williams  1 1971,  p.  14  5 J . 

Simple  trapezoidal  dams:  W.  D.  L.  Finn  [1967];  S.  P.  Neuman  and 
Witherspoon  [ 1970];  Parkin  [1971];  Pettibone  and  Kealy  [1971] 

(inhomogeneous  dam);  Fenton  [ 1972a]. 

Polygonal  dams:  W.  D.  L.  Finn  [1967]  (trapezoidal  dam  with  two  sheetpiles); 

R.  L.  Taylor  and  Brown  [ 196  7]  (with  rock  toe);  Volker  [1969] 

(with  cut-off  wall);  Neumann  and  Witherspoon  [ 1970]  (with  rock 
toe;  with  toe  drain);  Kealy  and  Busch  [1971],  Kealy  and  Williams 
[1971,  p.  152], 

Dams  with  general  geometry:  Nonlinear  flow  through  a dam  with  curved 
downstream  face,  and  nonlinear  flow  through  a polygonal  dam 
with  cut-off  wall  (Volker  [1969]). 

Circular  channel  in  a finite  layer  (pond):  Neuman  and  Witherspoon  [1970]. 

Single  fully-penetrating  well:  Neumann  and  Witherspoon  [1970]. 

Single  partially-penetrating  well:  R.  L.  Taylor  and  Brown  [1967]. 


2.3.  Trial  free  boundary  methods  with  Galerkin  methods 

By  a Galerkin  method  we  understand  any  method  which  involves 

(k) 

approximating  the  solution  u (x)  by  a known  smooth  function  F(c,x) 
depending  upon  an  N- vector  c which  is  determined  by  minimizing  the 


error  in  some  sense. 
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References  on  Galerkin  methods  include:  Kantorovich  and  Krylov  [ 1958]  , 
Cavendish,  Price,  and  Varga  | 1969]  , Mikhlin  [ 1971]  , Birkhoff  [ 1971]  , Varga 
[1971]  , 

In  all  the  applications  of  Galerkin  methods  to  trial  free  boundary  methods 

of  which  we  are  aware,  the  governing  differential  equation  has  been  Laplace's 

(k ) 

equation  and  the  approximate  solution  ul  has  been  represented  in  the  form 


UW(X,  = X,  . 


P„ti) 


N 

y 

LJ 

j = l 


cfL/x), 


(1) 


or  in  an  equivalent  form,  where: 

(i)  Each  function  p (x),  0 < j < N,  is  a solution  of  Laplace's 

equation. 

(ii)  F(c,  x)  satisfies  some  of  the  boundary  conditions  on  the  known  part 

of  for  all  c. 

(b) 

(iii)  The  function  p (x)  incorporates  the  expected  singularities  of  u '( x ). 


(k) 

The  values  of  the  coefficients  c.  have  been  determined  by 
satisfying  the  remaining  boundary  conditions  in  one  of  two  ways: 


(a) 


Collocation:  The  boundary  conditions  /9u 


(k) 


0 are  imposed 


at  N points  thus  leading  to  a system  of  N equations  for  the 

.(k) 


coefficients  c 


j • 
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(k) 

(b)  Least  squares:  The  boundary  conditions  0 are  imposed  at 

M > N points,  and  the  resulting  system  of  M equations  is  solved 
by  least  squares . 

We  now  summarize  the  various  applications  of  the  Galerkin  method. 
Garabedian  [19  56;  1956a,  p.  664]  considers  the  axially  symmetric 
Riabouchinsky  cavity  between  circular  disks.  Garabedian  takes  great 
care  in  the  choice  of  the  functions  p. . He  uses  least  squares  to 

/ Vr  \ 

determine  c^  with  N 10,  M - 24  (Garabedian  [ 1956a])  and  N 20, 
M 62  (Garabedian  [ 19  56]  ). 

Chappelear  [1961]  considers  periodic  progressing  gravity  waves  in 
water  of  uniform  depth.  Chappelear  takes 


ufk’(x,y)  -c^'x4  X c!klsin(2irjxA)cosh(2-njyA)/(2TTi/V) 

“h  ’ 0 ft  J 

where  \ is  the  wavelength.  Chappelear  uses  least  squares  to  determine 
Ik ) 

c and  it  appears  from  the  graphs  in  the  paper  that  M = 13. 

Kirkham  [ 1964]  considered  the  problem  of  axisymmetric  porous 

flow  towards  a fully-penetrating  well  of  radius  r in  a region  of  influence 

w 

of  depth  h and  radius  r . Kirkham  sets 
e e 

(Vy'^oMH  X cfk,c0sh(^]Co(Umr/rw'/c°sh(^]  > 

m=l  1 ' w w ' 


where 
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u i u r u r u r u r 

,,  / in  , , , m , ,,  . in  (•  , / m c\  ,,  , m , 

°o<r>  =Jo'7'>  V~ '"■>  • Vt-1  Vr1  ■ 

WWW  W V/ 

the-  functions  T and  Y arc  Bessel  functions,  and  the  constants  u are 
• 0 0 rn 

the  roots  of  C (u  r /r  I 0.  The  function  p (r,y)  is  a known  function 
0 m e w u 

expressed  as  an  infinite  series.  Kirkham  set  N = 4 and  used  collocation 
to  determine  c It  should  be  observed  that  Kirkham  explains  his  method 
in  terms  of  a "fictitious  region",  but  this  in  fact  is  merely  a means  of 


justifying  the  approximation  for  u^  . 

Mason  and  Farkas  [ 1971,  1972]  consider  the  problem  in  which 
fresh  water  is  supplied  to  a canal  in  the  center  of  a long  island  with  an 

impervious  surface.  Mason  and  Farkas  use  conformal  mapping  to  straighten 
out  1 a corner  in  the  original  domain  and  thereby  remove  a singularity  in 
the  solution.  In  the  transformed  z x + iy  domain  the  solution  is 
sought  in  the  form 


4>  '^(z)  + a >v  '(z)  = i + Yj  c\  z 

j=l  J 


(k ) j —1/2 


Collocation  is  used  to  determine  cv  . The  collocation  points  are  chose^i 
in  a special  way  which  is  explained  and  justified  by  Mason  [ 1969]  . Com- 
putations were  carried  out  for  N = 8,  10,  12,  and  14. 


We  conclude  with  two  remarks: 

(i)  It  is  not  necessary  for  the  function  F to  have  the  form  (1),  nor  is  the 
Galerkin  method  necessarily  restricted  to  Laplace's  equation. 
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(11)  The  Galcrkin  method  .has  several  advantages: 

(a)  It  is  easy  to  handle  the  trial  TBS  r 

(b)  The  approximate  solution  obtained  is  in  a form  which  facilitates 
error  estimates. 

(c)  The  value  of  Cu  is  readily  obtained. 

We,  therefore,  believe  that  the  Galerkin  method  deserves  greater  attention 
puiticularly  for  problems  with  Laplace's  equation  and  relatively  simple 
geometry  for  which  good  choices  of  F(c,  x)  can  be  made. 

2.4.  The  trial  free  boundary  method  with  integral  equations  and  surface 
singularities 

One  of  the  basic  methods  of  solving  a fixed  boundary  value  problem 

is  to  reformulate  the  problem  as  an  integral  equation.  It  is,  therefore,  very 

, (k) 

natural  to  use  this  approach  to  compute  u^  . 

The  first  use  of  integral  equations  in  connection  with  trial  free 

boundary  methods  is  due  to  Trefftz  who  in  his  thesis  of  1914  (the  greater 

part  of  which  appears  in  Trefftz  [ 1916]  ) considered  the  problem  of  an 

inviscid  axially  symmetric  wall  jet.  Since  then  integral  equations  have 

been  used  by  several  authors.  The  approach  has  not  been  as  widely  used 

as  might  be  expected  for  a number  of  reasons: 

(i)  In  general,  the  numerical  solution  of  the  integral  equation  requires 

the  solution  of  a non- sparse  system  of  linear  algebraic  equations,  and 
this  only  became  convenient  after  the  introduction  of  computers.  (Because 
of  the  special  features  of  his  problem  Trefftz  was  able  to  avoid  this.) 
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(ii 1 The  approach  is  most  useful  in  axially  symmetric  problems,  for 
which  complex  function  techniques  are  not  available.  In  axially 
symmetric  problems,  however,  the  kernels  of  the  integral  equations 
are  non-trivial  elliptic  integrals. 

(iii  A particular  FBP  can  be  reformulated  as  an  integral  equation  in 

several  ways.  Moreover,  the  integral  equations  are  often  derived 
in  rather  tortuous  ways  using  physical  ideas.  In  consequence, 
the  literature  is  very  fragmented. 

Despite  these  difficulties,  integral  equations  are  being  used  increasingly 
often,  and  this  trend  will  probably  continue. 


2.4,1.  The  method  of  Trefftz 

The  method  described  here  was  first  used  by  Trefftz  [1914,  1916]  to 
solve  the  problem  of  an  inviscid  axially  symmetric  wall  jet.  Gilbarg  [i960, 
p.  431]  gives  a useful  summary  of  the  work  of  Trefftz. 

The  geometry  of  an  axially  symmetric  wall  jet  is  shown  in  Figure  1. 
Applying  Green's  identity  we  find  that 


2it  cf  ( 0 ) = + J 

86. 


4>M 


_8 

8n 


r 

g pa 


dS 


J 

86. 


1 9e (q ) 

r 8n 
Pa  q 


dS 


where  6^  is  the  (three-dimensional)  axisymmetric  FD,  p and  q are 

points  on  86  , r is  the  distance  between  p and  q,  and  <j-  is  the 
j pq 


velocity  potential. 
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Figure  1 : An  axially  symmetric  wall  jet 


/.  ^ is  axially  symmetric.  Let  £ be  the  cross-section  of  in  the  xy- 
plane.  £ is  symmetric  with  respect  to  the  y-axis,  and  we  denote  by  3 £ the 
boundary  of  £ lying  to  the  right  of  the  y-axis.  Then  equation  (l)  can  be  reduced 
to  an  equation  in  the  xy- plane  of  the  form 

z-n- 4 ( s ) = I K.(s,t'<Htidt  + j K_(s,t1  dt  , (2' 

J.  “t  „ c on 

d£  3£ 

where  s and  t denote  arc-length  along  3)9,  and  and  are 

kernels  which  depend  upon  s,  t and  the  FB  r (the  boundary  of  the  jet l. 
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More  precisely,  if  the  coordinates  of  the  points  s and  t are 


(xis  ,y(s"  and  (x(t),y(t'l,  respectively,  then 


K,(s,tl  - 2x(ti  7 — V s,t)  , 
1 dn  ’ ’ 

t 


and 


K (s,ti  = -2x(t)V(s,t)  , 

where 


l 

V(s,t)  - - - J — 

0 {[  y fs)  - y (t) ] 2 + f x(sl]  2 + [ x(tl]  2 - 2x(s)x(t)  cos  oo}1^2 

To  avoid  dealing  with  an  infinite  integral  the  jet  is  truncated  at 

a line  y -x,  so  that  the  boundary  consists  of  three  components 

namely  the  1’B  1',  the  part  d $ of  the  line  y 0,  and  the  segment 

d of  the  line  y - (see  Figure  1). 

I'  and  a,  A are  streamlines  so  that  ^ 0 on  these  curves.  It 

1 an 

is  assumed  that  the  velocity  is  constant  on  a^jt,  and  the  problem  is  so 

normalized  that  ~ = 1 on  0 _ jS  . Thus,  for  given  F the  second 
an  2 

integral  in  equation  (2)  is  a known  function  of  s,  f(s;D  say, 
and  the  equation  takes  the  form 


2tt<j>(s ) = f(s ; F)  + f K (s,t)b(t)dt  . (31 

d&  1 

A further  simplification  arises  because  K^(s,t)  = 0 if  both  the  points 

s and  t lie  on  the  fixed  boundary;  this  is  best  seen  from  equation  (1) 

ar 

PQ 

by  noting  that  if  p and  q both  lie  on  the  plane  y - 0 then  - 0. 

9nq 

Thus, 
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K (s , t)p(t)dt,  for  sc  9^ 


2tT(p(s)  - f(s;l  ) + f 

Iua/ 

2tt<>(s)  f ( s ;1 ' ) 1 / K (s,t)*(t)dt,  for  st  r u 9^  . 

The  additional  condition  on  the  FB  is  that  o(s)  - aQ  + s on  I , 

where  o is  the  value  of  6 of  the  point  of  separation.  It  follows  that 

a a + f on  9„A  where  f is  the  length  of  1 . 
p <-*0  2 

The  procedure  follows  by  Trefftz  is  as  follows: 

(i)  Guess  the  FB  T ^ . 

(ii)  Compute  the  value  of  6 on  the  known  boundary  9^  by  means  of 


2u<,(1|<s)  , Ks-.r1'1)  * f K1(S,t)|4)H)dt+  / aKl(=,tKVnUt, 

rai  !,2” 


.(1) 


for  st  9^®  • 


(Hi)  Compute  the  value  of  £ on  1 by  n.eans  of 


2tt6(1)(s)  ^ f(s;l  (1))  - / K1(s,t)4>(1)(t)dt  + 

9 } S 


/ K^s,  t)(*0+f  )dt  + /(1)K1(s,t)[  «i,0  + tldt. 


for  s t I 


(1) 


9/ 


(iv)  If  (j)^^ (s ) - + s on  then  the  problem  is  solved.  Other- 

wise we  must  adjust  rQ).  Trefftz  [ 1916,  p.  47]  suggests  the 
following  strategy: 

(a)  if  6(1)(s)  > o0  + s on  r(1),  then  r(1)  should  be 

moved  away  from  the  axis; 

(b)  if  <|>^(s)  < 6Q  + s on  rfI\  then  1 (1)  should  be  moved 


towards  the  axis, 
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Trefftz  used  power  series  expansions  and  a planimeter  to 

approximate  the  integrals.  Because  of  the  effort  involved,  Trefftz  only 

carried  out  two  iterations  and  he  estimated  that  .60  < C < .62. 

— c — 

Steffen  [ 197  3]  has  implemented  a modification  of  Trefftz 's  method 
on  a computer.  Steffen  chooses  n values 

-k  < y <y  ,<•••<  y,  < 0 
n n-1  1 

and  defines  the  kth  approximate  FB  r^'  to  be  the  quadratic  (sic)  spline 
](£ 

x = x (y)  passing  through  the  points 

(1,0).  (x[,y,> (*[,yn>,  - -2H(-i  - -) 

(k) 

which  has  a horizontal  tangent  at  (1,  0).  Given  1 Steffen  computes 
6^  as  done  by  Trefftz  and  then  determines  the  error 


Fk(y)  = tj,(k'(xk(y),  y)  - (6Q  + fk(y)) 


k k 

where  f (y)  is  the  length  along  r expressed  as  a function  of  y. 

The  formulas  of  Steffen  [ 1973,  p.  14]  differ  from  those  of  Trefftz  in  two 
minor  respects:  the  formulas  are  expressed  in  terms  of  x and  y 
rather  than  s and  t;  and  Steffen  makes  use  of  the  identity  (see 
equation  (1)) 


2tt  = f (— — )dS 

J an  r 


To  find  the  n values  x 


scheme: 


k+1 
1 ’ 


q pq 

kfl 


Steffen  uses  the  iteration 
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M 


xk)1  = -c  V Fk(y, ) I (i  - t + 1)  rL  + xk  for  1 < j < n, 

1 t l ( 

where  c and  L are  constants  chosen  to  improve  convergence:  the 
best  values  were  found  to  be  c - .031  and  L = 1.  Steffen  [ 1973,  p.  16] 
reports  that  the  iterations  did  not  converge  completely,  but  oscillated. 
Steffen  [ 1973,  p.  16]  also  observes  that  he  tried  Newton's  method 
without  success. 

At  this  point  it  is  appropriate  to  emphasize  the  essential  features 
of  the  method  of  Trefftz: 

(i)  The  method  is  based  upon  equation  (1)  which  involves  both  <*> 

and  , and  use  is  made  of  the  fact  that  both  6 and 

3n  ’ On 

are  known  on  the  FB. 

(ii)  There  is  no  need  to  solve  an  integral  equation  approximately 
because  the  geometry  is  such  that  the  kernel  K^(s,  t)  vanishes 
on  the  fixed  boundary. 

For  many  years  the  work  of  Trefftz  was  the  basic  work  on  axially 
symmetric  problems.  The  method,  or  variations  thereof,  was  used  to 
solve  the  following  problems: 

(a)  Axially  symmetric  jet  impinging  normally  on  a plate:  Schach  [ 1935]. 

(b)  Wedge  penetrating  a fluid:  Wagner  [1932,  p.  213],  Pierson  [ 1950]. 
(As  Pierson  [19  50,  p.  2]  remarks,  the  work  of  Wagner  is  brief 

and  cryptic.  Pierson  [ 19  50,  Appendix  I]  gives  more  details 
but  is  also  not  readily  comprehensible.) 
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(c)  Cone  penetrating  a fluid:  Shiftman  and  Spencer  [ 19  51]  (the 

computations  were  performed  by  Hillman  in  1946  in  a report  that 
we  have  been  unable  to  obtainC 

(d'  Axisymmetric  magnetosphere  under  uniform  external  pressure: 

Slutz  [ 196  2]  . 

The  paper  of  Shiftman  and  Spencer  [19  51]  is  very  clear.  The 
method  used  differs  from  that  of  Trefftz  in  two  respects:  the  kernel 
Kj(s,t\  does  not  vanish  on  the  fixed  boundary  so  that  it  is  necessary  to 
solve  an  integral  equation  to  obtain  the  values  of  <f  on  the  fixed 
boundary:  the  FB  is  moved  by  a global  method. 

The  work  of  Slutz  [ 196  2]  is  an  interesting  contribution  to  a recent 
problem.  There  is  no  fixed  boundary  so  that  it  is  not  necessary  to  solve 
an  integral  equation.  The  FB  is  moved  by  a global  method. 

2.4.2.  The  trial  free  boundary  method  with  integral  equations 

A classical  method  for  solving  the  Dirichlet  problem  for  Laplace's 
equation 

u + u =0,  in  S 

xx  yy  (1) 

u - f,  on  96 

is  to  represent  the  solution  u as  the  potential  of  a double  layer  of 

density  p on  96;  that  is 

u(x)  = J F(x,|_(t))p(t)dt  , (2) 

96 
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where  t denotes  distance  along  <>ft,  the  density  n is  an  unknown 
function,  ^ £_(t'  is  the  point  on  3ft  corresponding  to  t,  and 


/ >. , 1 cos  <} 


where  <r  and  r are  as  in  Figure  1. 


For  fixed 


so  that  F(x,£_)  is  a solution  of  Laplace's  equation.  Thus  the  integral  in 


(2),  being  the  superposition  of  solutions,  is  also  a solution  of  Laplace's 


equation.  The  term  double  layer  refers  to  the  fact  that  F can  be  inter- 


preted as  the  potential  due  to  a double  layer. 


In  order  to  check  that  u satisfies  the  boundary  conditions  it  is 


necessary  to  study  the  behaviour  of  the  integral  in  (21  as  x approaches 


the  boundary.  If  the  boundary  is  smooth  at  the  point  t then  it  can  be 


shown  that 


lim  J F(x,£.(t))p(tldt  (j.(s ) + J K(s,t)p(t)dt 
_x-*£(s'  eft  “aft 


where 


K(s  t)  = - ~ arc  tan 

tt  dt  4 (t)  - £(s 


Thus,  u given  by  (2)  will  solve  the  boundary  value  problem  (1)  if  |jl 


satisfies  the  Fredholm  integral  equation  of  the  second  kind, 


M-(s)  + J K(s,t)p(t)dt  - f(s)  . 
3ft 


I 


W 


it  is  possible  to  rewrite  the  above  results  in  another  form.  The 
expression  (2)  is  equivalent  to  the  Stieltjes  integral 

u(x)  ~ J (jl ( t )dtu 
80 

where  co  is  shown  in  Figure  1 . The  limit  (5)  is  readily  seen  t o be  a 
conseouence  of  (8i.  Finally,  equation  (7)  is  equivalent  to  the 
Stieltjes  integral  equation 

pis)  + ^ J p(t )d  v(s,t)  f(s)  , 


(8) 


19) 


where  is  as  in  Figure  1. 

The  derivation  of  the  preceding  results  is  described  in  many  textbooks 
including:  Kantorovit  ch  and  Krylov  [ 19  58]  , Courant  and  Hilbert  (1962,  p.  298], 

The  integral  equation  (91  must  be  solved  by  some  numerical  method. 
General  references  on  the  numerical  solution  of  integral  equations  include: 

AnseF  ne  [ 197  1 ] . Noble  [ 1971]  . Delves  and  Walsh  [ 1974]  , Atkinson  [ 1976]  , 
Baker  [ to  appear]  . 

The  numerical  solution  of  equation  (91  is  considered  by  Kantorovich 
anc  Krylov  [ 19  58]  , Symm  [ 1964]  , Jaswon  and  Symm  (to  appear]  . Cryer 
[ 1970,  p.  103]  gives  further  references . 

There  is  a close  relationship  between  equation  (9)  and  various 
integral  equations  which  arise  when  constructing  conformal  maps:  see 
Gaier  [ 1964]  , Symm  [ 1966,  1967.  1969]  . 
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that  if  the  matrix  is  too  large  to  be  held  in  the  main  store,  then  the 
computations  must  be  arranged  so  as  to  minimize  transfers  between  the 
main  store  and  the  secondary  store;  an  efficient  way  of  doing  this  is 
described  in  a seldom -quoted  paper  of  Barron  and  Swinnerton-Dyer  [ 1960/61] 
which  was  intended  for  a magnetic-tape  secondary  store  but  is  applicable 
to  other  types  of  secondary  store.  The  time  estimates  given  by  Barron 
and  Swinnerton-Dyer  [ 1960/61]  are  too  low  because  only  the  time  required 
for  data  transfer  is  considered.  When  allowance  is  made  for  the  computation 
time,  the  time  estimates  are  very  accurate  (Cryer  [ 1962,  Chapter  5,  p.  7]). 

Since  the  system  (10)  is  not  sparse  it  might  be  thought  that  the 
solution  of  (10)  would  require  much  more  computation  time  than  the  solution 
of  the  sparse  linear  equations  arising  when  finite  elements  or  finite 
differences  are  used  to  solve  (1).  However,  equation  (9)  is  a one- 
dimensional equation  while  equation  (1)  is  a two-dimensional  equation. 

Thus,  for  comparable  accuracy,  the  order  of  the  system  (10)  is  far  smaller 
than  the  order  of  the  system  arising  from  finite  differences  or  finite  elements. 
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We  now  indicate  some  related  results: 

(i)  Three  dimensions 

The  results  can  be  extended  to  any  number  of  dimensions.  In 
three  dimensions 


F<X,ii  ~ T3 

’ / TT  1 1 


1 


2ir  6n  |x  - 4.  | 

and  the  kernel  K(s,t)  is  modified  appropriately.  See  O.  D.  Kellogg  [ 19  53, 
p.  28o]  , Garabedian  [ 1964 , p.  334]. 

General  programs  for  solving  the  resulting  integral  equations  have 
been  developed  by  Hess  and  Smith  [ 1967]  and  Hess  [ 1973]  and  applied  ' 
to  the  computation  of  the  flow  past  obstacles  such  as  aeroplanes  and 
mis  siles . 

For  axisymmetric  problems  the  axisymmetry  can  be  used  to  reduce  the 
integral  equation  to  a one-dimensional  equation.  The  surface  singularities 
are  often  called  ring  sources  or  ring  vortices.  Unfortunately,  the  kernel 
of  the  one-dimensional  equation  involves  elliptic  integrals  and  this 
increases  the  computation  time.  See  also  Riegels  [ 19  52]  , Landweber  [19  51, 
19  59]  . 

(ii)  Hon- smooth  oft 

If  8ft  has  an  exterior  corner  of  angle  a(s)  at  the  point  s on  o& 
then  (5)  must  be  modified: 


lim  J F(x,£(t))p(t)dt 
x--^(s)  8* 

The  integral  equation  (71  becomes 


Q(S) 

TT 


M-(s)  + J K(s , tlpltldt 
3ft 


“ p(s)  + J K(s,ty(t)dt  f ( s ) . 
8fl 
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However,  and  this  is  its  great  advantage,  equation  (9)  remains  valid, 
that  is 

fj.(s ) + ^ J ji(t)d^(s,  t)  = f(s)  . (11) 

Equation  (11)  was  first  considered  by  Carleman,  but  Radon  [1919] 
developed  a much  more  satisfactory  theory.  Recently  the  work  of  Radon 
has  been  extended  by  Burago,  Krai,  and  others:  Radon  defined  eft  to  be  of 
bounded  rotation  if  the  slope  0(s)  (see  Figure  1)  is  of  bounded  variation. 
Radon  showed  that  if  8ft  is  of  bounded  rotation  and  has  no  cusps  then 
(9)  has  a unique  continuous  solution  p for  each  continuous  f.  Cryer  [ 1970] 
gives  a detailed  discussion  of  the  theory  and  a comprehensive  bibliography. 

The  general  theory  of  Radon  is  not  necessary  in  the  present 
context,  but  most  FBPS  do  involve  boundaries  Eft  with  corners,  so  that 
the  theory  which  is  usually  given  in  textbooks  which  assumes  a smooth  8ft 
and  is  based  on  equation  (7)  is  not  adequate.  Equation  (11)  has  been  used 
with  trial  free  boundary  methods  by  Cryer  [ 1962,  Chapter  1,  p.  11]  , and 
Symm  [ 197  5]  . Slutz  [ 1962]  used  the  same  approach  in  a three-dimensional 
magnetohydrostatic  problem  which  he  solved  by  the  method  of  Trefftz  (see 
section  2.4.11;  an  interesting  aspect  of  the  problem  treated  by  Slutz 
is  that  the  FB  has  a corner. 

(iii)  Single-layer  potentials 

Instead  of  representing  u(x)  as  the  potential  of  a double  layer 
of  density  p,  one  can  represent  p as  the  potential  of  a single-layer 
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of  density  p.  For  two-dimensional  problems  we  have 

u(x)  J p(t)(n|x  - t > | dt  , (12) 

b& 

instead  of  (4) . 

If  the  boundary  condition  is  a Neuman  condition,  then  equation  (12) 
leads  to  a second  order  Fredholm  integral  equation  for  p with  the 
kernel  K(s,t)  given  by  (6).  However,  if  the  boundary  condition  is  a 
Dirichlet  condition,  then  (12)  leads  to  the  first  order  equation 

f ( s ) J p(t)ln[|(s)  - |(t)  | dt  . (13) 

Conventional  wisdom  has  it  that  equations  such  as  (13)  should  be 
avoided.  However,  it  has  been  found  by  Symm  [1964]  , Jaswon  and  Symm 
[to  appear]  , and  others  that  equation  (13)  can  be  solved  numerically  without 
difficulty . 

So  far  as  we  are  aware,  a rigorous  treatment  of  the  numerical 
solution  of  singular  equations  of  the  form  (13)  is  not  available.  There  is 
a considerable  literature  on  the  numerical  solution  of  equations  of  the 
second  kind  with  singular  kernels  of  which  we  merely  mention  Hertling  [ 1971] 
and  Mclnnes  [ 197  2]  (who  generalizes  the  theory  of  collectively  compact 
operators  to  the  case  of  linear  singular  integral  equations  in  a Holder 
space.  ) There  is  also  a considerable  literature  on  the  solution  of  first 
order  integral  equations;  a recent  reference  is  Strand  [ 1974]  . 
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1 1 v C in.binc  : potentials 


It 


is  possible  to  represent  u(xi  in  the  form 
u(x‘  ~ J n(t)  — ( n |x  - £(t)  Idt  + J n(t)fn 


- |(t)  Idt  , 


where  U r ft  ffl.  We  understand  that  this  approach  was  used  by 

Symm  [ 197  5]  to  solve  the  problem  of  porous  flow  through  an  earth  dam. 

(v ' Biharmonic  equations 

Symm  [ 1964]  and  laswon  and  Symm  [to  appear]  show  how  the 
biharmonic  equation  can  be  solved  using  integral  equations.  This  could 
be  an  effective  method  for  solving  FBPS  for  the  biharmonic  equation. 

Cvi)  Connections  with  complex  function  theory 

The  solution  u of  (1)  is  the  real  part  of  an  analytic  function 
w w(zi  - w(x  + iy).  There  is  a very  close  connection  between  the 
representation  (2)  and  the  representation  of  w as  a Cauchy  integral 


w(z ) 


1 

2TTi 


a® 


Mil 

i - z 


de . 


The  behavior  of  integrals  of  the  form  (14)  is  extensively  studied  in  the 
literature  on  complex  function  theory.  See  Muskhelishvili  [ 19  53a]  and 
Goluzin  [1969]  . For  applications  to  the  theory  of  elasticity  see 
Muskhelishvili  [ 19  53]  . 

The  use  of  integral  equations  for  trial  free  boundary  methods  has 


several  advantages: 


I 


\ JS.  I 

13'  General  domains  f,  can  be  easily  handled. 

(k  i 

ib  The  appr  ximate  solution  u,  is  expressed  in  terms  of  an  integral 

h 

( ^ ) 

involving  an  approximate  density  fj.,  so  that  the  values  of  u^ 

and  its  ierivatives  can  be  readily  computed  at  any  point  of 

( k 1 ( | 

Thus  the  boundary  condition  Cu  ‘ 0 on  d£  ' is  easily  checked, 

ir  C rner  singularities  are  easily  handled  by  increasing  the  density  of 
the  meshpoints  near  the  corners.  For  another  approach  see 
Symm  [ 197  3]  . 

(d  i Finite  difference  methods  and  finite  element  methods  immediately 

give  values  of  the  solution  at  mesh  points  , but  if  integral  eguations 

(ki 

are  used  then  the  values  of  the  solution  at  points  in  & must  be 
computed,  which  is  a straightforward  but  time-consuming  computation 
which  can  require  as  much  time  as  the  solution  of  the  integral 
equation.  For  fixed  boundary  value  problems  this  additional  computa- 
tion is  a definite  disadvantage  of  the  method  of  integral  equations. 
However,  for  trial  free  boundary  methods  the  values  of  the  solution 
in  ft  are  only  needed  (if  at  all)  after  the  final  iteration,  so  that 
the  additional  computation  is  only  a slight  disadvantage. 

The  main  disadvantage  associated  with  the  method  of  integral 
equations  is  that  the  method  is  limited  to  a few  special  differential 
equations  such  as  Laplace's  equation.  This  is  a severe  restriction,  but 
it  is  not  as  severe  as  might  be  thought:  we  have  compiled  an  extensive 
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bibliography  of  FBPS,  and  a very  large  number  of  these  FBI'S  involve 


Laplace's  equation. 

The  trial  free  boundary  method  with  integral  equations  has  been  used 
to  solve  the  FBPS  listed  below.  Information  about  the  type  of  surface 
singularity,  the  type  of  integral  equation(s),  and  the  numerical  method  used 
to  solve  the  integral  equation(s)  is  given  in  brackets. 

Axisymmetric  Riabouchinsky  cavity  flow:  Armstrong  and  Dunham  ( 19  53] 

(ring  vortices). 

Axisymmetric  cavity  flow  (Helmholtz  model  for  cones  and  spheres: 
Riabouchinsky  model  for  disks):  Struck  [ 1970]  (ring  sources; 
equations  of  first  and  second  kind;  Nystrom  method). 

Axisymmetric  wall  jet:  B.  W.  Hunt  [ 1967,  1968]  (ring  vortices;  equation  of 
second  kind;  Krylov  and  Bogoliubov  method). 

Porous  flow  through  rectangular  dams:  Symm  [ 197  5]  (single-layer 

potential'.  Cryer[l962,  chapter  1,  p.  11]  used  a double-layer 
potential  to  check  the  solution  of  Shaw  and  Southwell  [1941]  but 
did  not  iterate. 

Elastic-plastic  torsion:  Ponter  [ 1966]  and  Jaswon  and  Ponter  [ 1963] 

(double  layer;  equation  of  second  kind). 

2.4.3.  The  self-consistent  method  and  the  moment  method 

In  this  section  we  describe  two  methods  which  use  surface  singularities 
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The  self-consistent  n.Gthod  was  introduced  by  Mead  and  Beard  [ iv64j 


to  compute  the  earth's  magnetopause  for  the  symmetric  case,  and  was 
later  used  by  Olson  [ 1968,  1969]  for  the  general  case. 

The  problem  is  as  follows.  A stream  of  charged  particles,  called 
the  solar  wind,  flows  from  the  sun  towards  the  earth.  The  solar  wind  has 

A.  N 

density  p and  velocity  Vv  where  v is  a unit  vector.  The  eart  . is 
represented  as  a magnetic  dipole  of  strength  M.  The  meridian  plane 
is  the  plane  containing  the  earth's  dipole  and  v.  The  equatorial  plane 
contains  v and  is  perpendicular  to  the  meridian  plane.  The  earth's 
magnetic  field  deflects  the  solar  wind  and  creates  a three-dimensional 
region  £ near  the  earth  which  is  not  penetrated  by  the  solar  wind.  The 
two-dimensional  boundary  r of  fl  is  a FB,  the  magnetopause . 

The  governing  equations  can  be  set  up  by  noting  that: 

(i)  There  is  a magnetic  field  inside  the  magnetopause  I which 
satisfies  Maxwell's  equations. 

(ii)  There  is  no  magnetic  field  outside  r. 

(iiit  The  solar  wind  particles  must  behave  in  a prescribed  fashion  on 

r.  The  most  common  assumption  is  that  there  is  specular  reflection; 
that  is,  that  the  particles  are  reflected  by  the  magnetopause 
in  the  same  way  that  light  would  be  reflected  by  a mirror  r. 

In  terms  of  appropriate  polar  coordinates  (r,  0,  c the  FB  is 
defined  as  the  surface 

r R(0,c)  . 
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The  equatorial  cross-section  of  the  FB  is  smooth  while  the  meridional 


cross-section  contains  two  discontinuities  - the  neutral  points  where 
the  magnetic  field  is  zero. 

The  total  magnetic  field  inside  1 may  be  thought  of  as 
consisting  of  two  components  namely  the  geomagnetic  dipole  field  B 

~g 

and  a field  B caused  by  surface  singularities,  namely  currents  flowing 
— c 

in  I'.  The  field  B is  known,  and 

~g 

B B r •)  B 0 4 B 6 
— c cr  c0  c'b 

where  r,  0,  and  <*>  denote  the  unit  vectors  in  the  coordinate  directions. 
11  specular  reflection  occurs  then  on  the  FB  we  have  the  condition 

(Olson  [ 1968,  p.  31]  ), 


Cu:  l(r-^S)0 


1 o R A A 

„ . . (—  )b)  X ( (B  -2  sin  0 sin  <b)r  + 

R DO  R sin  0 06  cr 


4 (B  4 cos  0 sin  6)0  + (B  4-  cos  6)6  I 
cO  C6 


- 8 it p VM  1R3(cos  0 + (~^))  , 


(1) 


= 0 . 


The  current  field  B satisfies  the  equation  (Olson  [ 3968,  p.  35]), 
- c 


. (n  XB(q))Xr 


(2) 


(r  ) 

-pq 


where  n is  the  unit  normal  to  r at  the  point  q.  where  dS  is  a 

q q 

surface  element  at  q,  and  where  r is  the  distance  from  p to  q. 

~pq 
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In  the  self-consistent  method  the  field  in  guessed  -■  usually 

R*0>  B . Equation  (1)  is  then  integrated  using  finite  differences  to 
-c  -g 

determine  I'<0'.  Next,  equation  (2)  is  used,  with  Bg  + Bf>  , to 

determine  B^'  and  the  process  is  repeated.  The  self-consistent 
~c 

(k) 

method  is  thus  a trial  free  boundary  method,  the  trial  FBS  I being 
found  by  an  integral  method  (see  section  3.2.2).  The  governing  equation 
is  Laplace’s  equation  expressed  in  the  form  (2)  using  double-layer  surface 
singularities . 

In  addition  to  the  work  mentioned  above,  Baker,  Beard,  and  Young 
| lg64]  have  applied  the  self-consistent  method  to  some  magnetopause 
problems  with  known  solutions  and  Board  f 1964]  gives  a good  description 
of  the  method. 

So  far  as  we  are  aware,  the  self-consistent  method  has  only 
been  used  to  compute  the  magnetopause.  The  method  would  appear  to 
have  applications  to  many  other  FBPS,  particularly  problems  of  magnetic 
containment. 

The  moment  method  was  introduced  by  Midgley  and  Davis  [ 1962] 
to  solve  a magnetopause  problem  for  the  case  of  constant  pressure,  and 
was  subsequently  used  by  Midgley  ( 1963]  and  Midgley  and  Davis  [ 1963] 
for  the  case  of  specular  reflection.  The  moment  method  involves  the 
following  steps: 
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(1) 


The  FB  is  represented  in  parametric  form  in  terms  of  an  unknown 
vector  c (el.  Midgley  and  Davis  [1962]  assume  that  the  FB  is 
of  the  form 


while  Midgley  and  Davis  [ 1963]  assume  that  the  I'B  is  given  in  terms 
of  the  "flux  function" 

f(p  , <fc)  a(v)sin  tji  ’ V'(l  - [ ( 1 - u ) 2 + h(u,  v) ] 

where  v cos2<}>,  a(v)  is  half  the  radius  of  the  FB  at  z = -«, 
u p/n(v),  and  g and  h are  double  power  series  in  u and  v 
with  unknown  coefficients  c ■ . 

(ii)  Tot  an  assumed  FB  tire  surface  current  £_  is  determined  and 
the  resulting  vector  potential  A outside  the  magnetopause  is 
expanded  in  a series  of  functions  D which  are  related  to  the 
Legendre  polynomials, 

A ( p ) = f f f^r  hSa  = V I (c)D  (2). 

— J J [f  | q a — a 

I pq  « 

The  coefficients  I are  called  the  moments. 
a 

(iii)  The  condition  that  the  magnetic  field  outside  I'  should  vanish 
is  expressed  in  terms  of  the  moments  I . The  resulting  system 
of  equation  for  the  unknown  coefficients  c is  solved  iteratively . 
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The  moment  method  is  thus  a trial  free  boundary  method,  the  trial  tree 


boundaries  being  found  by  a global  method.  (See  section  3.  2.  3). 

Although  Midgley  and  Davis  did  obtain  results  they  experienced  consider- 
able difficulties  and  the  moment  method,  in  this  form  at  least,  does 
not  seem  to  be  a competitive  method. 

2.5.  Trial  free  boundary  methods  with  analog  methods 

Analog  methods  have  long  been  used  to  solve  1'BI’S  by  the  trial 
free  boundary  method.  Bear  | 197  2,  chapter  11]  gives  an  excellent 
discussion  of  analog  methods  for  porous  flow  problems,  and  Robertson 
[ 1965]  describes  analog  methods  for  fluid  mechanics  problems.  Basic 
texts  on  analog  computers  include  Karplus  [ 19  58]  , Volynskii  and 
Buhman  [ 196  5]  . 

In  the  discussion  below  we  do  not  consider  analog  methods  which 
merely  involve  scaling,  such  as  the  modelling  of  porous  flow  problems 
using  a "sand-box". 

Two  essential  differences  between  analog  and  digital  computers 

Accuracy.  In  an  analog  computer  it  is  difficult  and  expensive 
to  achieve  great  accuracy,  and  one  must  usually  be  content  with 
1%  accuracy.  In  a digital  computer,  far  greater  accuracy  is 
possible.  It  is,  however,  important  to  note  that  the  great  accuracy 
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are: 


(i) 


of  a digital  computer  may  not  be  fully  utilized:  very  few  FBPS 


have  been  solved  with  an  accuracy  greater  than  1%  , and  in  many 
practical  problems  the  data  is  not  known  to  1%  accuracy. 

(n)  An  analog  computer  is  special-purpose  while  a digital  computer 
is  general  - purpose . An  analog  computer  can  thus  sometimes  be 

much  faster  than  a digital  computer  on  the  class  of  problems 
for  which  the  analog  computer  has  been  designed. 

At  present,  analog  computers  have  a dowdy  and  old-fashioned 

image  compared  to  the  glossy  and  modern  imago  of  digital  computers. 

We  suspect,  however,  that  the  value  of  analog  computers  as  efficient 

special-purpose  devices  will  be  appreciated  more  in  the  future. 
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2.5.1.  Trial  free  boundary  methods  with  electrolytic  tanks 

An  electrolytic  tan!-,  is  a tank  filled  with  a conducting  medium  - 
the  electrolyte  - such  as  a copper  sulfate  solution.  The  electric- 
potential  in  electrolyte  satisfies  Laplace's  equation.  When  solving  a 
FBP,  the  FB  is  represented  by  an  adjustable  surface,  often  made  of  wax. 

It  is  possible  to  solve  two-dimensional  problems  involving  equations  of 
the  form 

i ,pu)  + w(pl" = 0; 

for  such  problems  the  tank  has  a movable  bottom  which  is  adjusted  so 
that  the  depth  h of  the  electrolytic  is  proportional  to  p. 

The  following  FBPS  have  been  solved  using  the  trial  free  boundary 
method  with  an  electrolytic  tank: 

Fluid  mechanics  FBPS 

Axially  symmetric  Riabouchinsky  cavity  behind  a disk:  Young,  Gates 

Arms,  and  Eliezer  [ 1955,  Figure  10]  and  Birkhoff  and  Zarantonello 
[ 1957,  p.  232]  give  some  details  of  a solution  obtained  by 
P.  Marchet  and  M.  Malavard  in  1949. 

Axially  symmetric  jet  from  a pipe:  Abul-Fetouh  [ 1949]  and  Rouse  and 
Abul-Fetouh  [ 19  50]  . 
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Axially  symmetric  jet  sti iking  a plate:  P.  G.  Hubbard  [1949]  and 
Leclerc  | 19  90]  . 

Weirs:  Hay  and  Markland  [ 1958]  obtained  results  for  the  flow  over  a 
weir  for  a large  number  of  geometries. 


Porous  How  TBPS 

Until  recently,  electrolytic  tanks  were  widely  used  for  studying 
two  dimensional  and  three  dimensional  porous  flow  problems.  Bear  [1972 
p.  706]  and  Polubaiinova-Kochina  [1962,  p.  463]  discuss  the  problem 
of  seepage  through  earth  dams. 

The  adjustment  of  the  FB  in  an  electrolytic  tank  is  quite 
laborious:  Rouse  and  Abul-Fetouh  [ 19  50  , p.  423]  remark  that  their 
solution  required  two  hundred  hours  of  work.  As  regards  the  accuracy 
attainable  with  an  electrolytic  tank.  Hay  and  Markland  [ 19  58,  p.  68  and 
p.  69]  state  that  their  results  are  correct  to  1%  , and  this  figure  seems 


to  be  generally  accepted. 


Trial  free  boundary  methods  with  the  Hele-Shaw  analog 


It  was  observed  by  Ilcle-Shnw  in  1898  that  if  a viscous  fluid 
flows  slowly  between  two  parallel  vertical  plates  which  arc  a distance 
b apart,  where  b is  small  then  the  constitutive  equation  reduces  to 


2 

V - (u,  v)  = - — (yb  /p)grad(y  + p/y), 

where  u and  v denote  the  average  velocity  of  the  viscous  fluid  in 
the  x and  y directions.  There  is  thus  a direct  analogy  between 
viscous  flow  between  parallel  plates  and  porous  flow.  This  is  of  course 
not  surprising  since  all  the  theoretical  derivations  of  Darcy's  law  in 
porous  flow  use  models  involving  viscous  flow  through  narrow  passages- 
Bear  [ 1972,  p.  687]  gives  a detailed  discussion  of  the  Hele-Shaw  analog 
which  has  been  applied  to  many  porous  flow  FBPS.  The  apparatus  for 
observing  viscous  flow  between  parallel  plates  is  often  called  a 
Hele-Shaw  cell . 

The  Hele-Shaw  analog  is  also  of  interest  because  it  provides  a 
link  between  porous  flow  FBPS  and  viscous  fluid  mechanics  FBPS  such 
as  the  two-dimensional  flow  of  a viscous  bubble  (see  Saffman  and 
Taylor  [ 1958]  ) . 
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2.5.3.  Trial  free  boundary  methods  with  resistance  networks 

In  sections  2.1  and  2.  2 we  considered  the  methods  in  which  the 
fixed  boundary  value  problem 


7u 


(kl 


0,  in  i 


(k, 


„ (kl  . _(k) 

Qu  0,  on  c&  , 

was  approximated  by  finite  differences  or  finite  elements,  the  resulting 
system  of  algebraic  equations  being  solved  by  hand  or  on  a computer. 

An  alternative  approach  is  to  replace  the  algebraic  equations  by  an 
equivalent  electrical  network  of  resistances. 

Resistance  networks  have  been  widely  used  for  porous  flow  FBPS 
(Bouwer  [ 1967]  , Bear  [ 1972,  p.  710]  ).  An  interesting  aspect  of  such 
applications  is  that  it  has  proved  possible  to  automate  the  movement  of 
the  trial  FBS: 

(i)  Karplus  [ 19  56]  automatically  solves  the  problem  of  a water-cone  by 
connecting  a series  of  DC  analog  computing  units  to  the  resistance 
network  along  the  network  boundary  corresponding  to  the  FB. 

(ii)  Herbert  and  Rushton  [ 1966,  p.  72]  automatically  solve  the  problem 

of  flow  through  a porous  dam  by  adding  transistors  which  disconnect 

(k) 

portions  of  the  network  where  the  condition  u > y is  violated. 
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2.5.4.  Trial  free  boundary  methods  with  conducting  paper  and  carbon 

Wyckoff  and  Reed  [ 1 3 5]  pioneered  the  use  of  conducting  paper 
to  solve  porous  flow  1 BPS.  They  used  graphite -coated  paper  to  solve 
the  problem  or  seepage  through  an  earth  dam  for  four  different  geometries: 
rectangular  dam  with  or  without  down-stream  water,  symmetric  trapezoidal 
dam  with  faces  of  slope  30°  and  4 5°.  In  the  approach  used  by  Wyckoff 
and  Reed  the  velocity  potential  0 in  the  porous  flow  is  represented  by 
the  electric  potential  in  a sheet  of  resistance  paper.  The  procedure  used 
is  shown  schematically  in  Figure  1 for  the  case  of  flow  through  a rectangular 
dam... 


Figure  1 : Resistance  paper  analog  of  porous  flow  through  an 

earth  dam  (based  on  Wyckoff  and  Reed  [ 1935.  p.  396]) 
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ln  Figure  J the  rectangul.it  sheet  of  graphite  - coated  paper  AEDF  is 
provided  with  terminals  which  maintain  parts  of  its  boundary  at  prescribed 
electric  potentials.  The  terminals  DLi  and  AB  are  good  conductors  and 
are  kept  at  potentials  e^  and  e , respectively,'  so  as  to  represent  the 
upstream  and  downstream  waters.  The  terminal  BF  is  a resistance  strip 
which  represents  the  seepage  face:  the  endpoints  T and  B are  kept  at 
potentials  e^  and  e^  and  the  potential  varies  linearly  between  F and 
B.  The  FB  is  the  curve  CD  and  is  obtained  by  cutting  the  sheet  of 
paper.  All  the  boundary  conditions  are  satisfied  except  the  second  boundary 
condition  on  the  FB.  namely  that  the  potential  should  be  linear  along  the 
FB,  and  the  FB  is  gradually  cut  away  until  this  boundary  condition  is 
satisfied. 

The  paper  of  Wyckoff  and  Reed  is  very  carefully  written.  They 
discuss  the  various  sources  of  error,  and  comparison  between  their  analog 
results  and  exact  analytical  results  shows  that  they  can  achieve  an 
accuracy  of  about  1% . 

Subsequently,  Childs  used  the  method  of  Wyckoff  and  Reed  to 
solve  several  porous  flow  problems:  drainage  of  rain  water  towards  a 
row  of  drains  (Childs  [19  43,  19 4 5 , 19 4 5a ] ) ; flow  down  a slope  interrupted 
by  channel  (Childs  [1946]);  the  Ghyben-Herzberg  lens  (Childs  [ 19  50]). 

I 
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The  use  of  graphite  paper  is  limited  to  plane  i'BPS.  Babbitt  and 
Caldwell  | 19  48]  studied  the  axisymmetric  problem  of  a fully -penetrating 
well  by  cutting  a thin  carbon  wedge,  applying  the  appropriate  electrical 
analogues  of  the  boundary  conditions,  and  then  cutting  the  upper  surface 
of  the  wedge  (which  corresponded  to  the  FBi  until  both  boundary  conditions 
were  satisfied. 

We  conclude  with  two  remarks: 

(i  1 In  an  era  of  high-speed  computers  the  use  of  an  analog  method 
based  on  graphite-coated  paper  may  seem  prehistoric.  However,  it  is 
not  unknown  for  large  amounts  of  computer  time  to  be  spent  to  obtain  results 
which  are  little  better  than  those  of  Wyckoff  and  Reed. 

(ii)  The  work  of  Wyckoff  and  Reed  [ 19  35]  and  Babbitt  and  Caldwell 
[ 1948]  provides  valuable  information  about  trial  free  boundary  methods  for 
porous  flow  problems:  it  is  clear  that  such  methods  must  be  very  stable: 
and  the  results  of  Wyckoff  and  Reed  suggest  that  if  the  initial  FB  is  too 
high  the  successive  trial  FBS  form  a monotonically  decreasing  sequence. 

2.6.  Trial  free  boundary  methods  with  graphical  methods 

The  earliest  approximation  methods  for  elliptic  equations  were 
graphical  methods,  and,  according  to  Runge  and  Willers  [1915,  p.  165]  , 
Forchheimer,  Blasius  and  von  Mises  all  considered  using  trial  free  boundary 
methods  in  conjunction  with  graphical  methods;  we  consulted  these  early 
references  but  found  that  none  of  the  workers  had  implemented  these  ideas, 
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which  was  of  course  quite  understandable  because  of  the  labor  required. 

1 "•  C ascigrande  [ 19 32,  1934]  and  A.  Casagrande  [ 1940]  , who  solved 
several  porous  flow  TBPS  for  dams,  were  apparently  the  first  to  solve 
i B!  S using  graphical  methods.  To  our  own  knowledge,  graphical  methods 
were  extensively  used  to  solve  porous  flow  problems  up  to  i960,  but  we 
are  not  aware  of  any  applications  to  FBPS. 

Shaw  and  Southwell  [ 1941]  were  aware  of  the  work  of  A.  Casagrande, 
so  that  there  has  been  a continuous  development  from  graphical  methods 
to  finite  difference  methods  and  then  to  finite  element  methods. 

m 


3 . Step  2:  Movement  of  the  boundary 

We  recall  the  definitions  of  Steps  1 and  2: 

( k ) ( k ) 

Step  1.  Given  I ' let  & ' be  the  corresponding  domain.  Compute 

(k  I , , . (ki  ,, 

an  approximation,  u,  say,  to  the  solution  u of  the  problem 

h 

*u<k>  = 0,  in  i,ki  , 

3uW  . 0,  on  . 

ik)  (k)  (k+1 i 

Step  2.  Given  I ' and  u ' compute  a new  trial  FB  F by 

(ki 

requiring  that  Cu^'  should  be  approximately  equal  to  zero  on 
T i , e . move  the  boundary  from  F to  1 

Various  aspects  of  Step  3 are  considered  in  the  following  subsections: 


choice  of  boundary  conditions:  movement  strategy;  a simple  convergence 
proof;  and  numerical  experience. 


3.1.  Choice  of  boundary  conditions 


In  most  FB  PS  the  bound. my  conditions  B and  C arise  in  a natural 


way  and  it  is  therefore  natural  to  use  them  without  modification.  It  is  the 
purpose  cf  this  section  to  point  out  that  it  may  be  helpful  to  modify  the 
boundary  conditions  on  the  FB  before  using  a trial  free  boundary  method. 


Consider  the  case  of  a one-domain  FBP  when  the  governing  equnt 
is  a general  second  order  elliptic  equation 


c?u  - a..u  + 2a._(u  + a n + b.u  + b_u  = 0 (1) 

11  xx  12  xy  22  yy  lx  2 y ' ' 


and  when  the  boundary  conditions  on  the  f'B  I take  the  form: 


Su 


P — 
1 3n 


C u = 


P, 


3u 


2 3n 


+ OjU  - Vj  = °* 

+ «2u  - Y2  = 0 , 


(2) 


where  p^,  «■. , and  y.  are  smooth  function  of  x and  y.  It  is  assumed  that 


h “l 


P2  tt2 


Y 0,  on  r,  (3) 


since  otherwise  the  boundary  conditions  on  r would  he  dependent  or 
contradictory. 

It  is  implied  by  (2)  that  the  boundary  condition 


9_u 

3n 


is  used  to  compute 


H a u - Y,  = 0,  on  r, 

1 1 

and  that  the  boundary  condition 


(*») 


P1  f)n  + V " *i  = °*  on  F> 

Jk) 

is  used  to  move  r 

The  role  of  (4)  and  (5)  could  of  course  be  interchanged.  There  are, 
however,  many  other  possibilities.  If  e_(x,  y)  are  smooth  functions 
satisfying 


S11  e!2 


°21  S22 


^ 0,  on  T,  (6) 
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then  instead  of  using  B and  r one  could  use  t .o  equivalent  operator:; 


V + e!2C>  °n  F’ 


6 = 


B 


on  ai-r, 


(7) 


C = e21S  + °22C  ’ °n  r- 


The  question  which  arises  is  whether  6 and  C can  be  chosen  so  as 
improve  the  convergence  of  a trial  free  boundary  method. 

Noting  (3)  we  see  that  two  possible  choices  for  8 and  C tr 


form: 


and 


8 | : — + v = 0, 

1 'F  an  Y11 


cr  U + Yiz  = 0> 


Bz  T : U + V21  = 0, 

c : ^ + 7 = 0 . 

^2  an  ' 22 


(«] 


(9) 


Of  these  two  possibilities  the  first  has  the  following  advantages: 
(i)  The  boundary  value  problem 


<7u,k)  = 0,  in  *,kl  , 


- (k)  _ (k  > 

e u = o , on  a&  , 


to 


ke  the 
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can  often  be  reformulated  as  a variational  problem  in  which 


the  boundary  condition 


9u^'1  n I (k) 

~ 4 'll  on  1 - 


appears  as  the  "natural"  boundary  condition. 

(ii)  When  u 1 has  been  computed,  both  u^  and  -p-  uP^  = -yn 
h ’ h 9n  h 11 

(k) 

are  known  on  r without  further  computations.  This  makes 
it  easy  to  use  linear  interpolation  to  find  the  next  trial 
TB  r(k+1)  on  which 


Cfu(k)  = uhk)  + V12  = °-  (10) 

A A 

In  contract,  if  the  operators  {/i  , d } wore  used  then  only 

tL  c. 

u[*  * is  known  on  1'^.  To  satisfy 


- /ki  „ (k) 

C2U  n‘TL~  + V?2  = 0. 
2 0n  22 


on  r 


(k+1) 


(lit 


using  linear  interpolation  it  would  be  necessary  to  compute 
, _JL  (k)  dZ  (k)  r(k) 

both  — u,  and  — - u,  on  i by  numerical  differentiation, 
an  h . 2 h 

Bn 

(iii)  In  general,  the  solution  of  a boundary  value  problem  involving 
Neumann  boundary  conditions  is  smoother  than  the  solution  of 
the  corresponding  problem  involving  Dirichlet  boundary'  conditions. 
One  might,  therefore,  hope  to  obtain  somewhat  greater  accuracy 

N 

A A A A 

using  c^}  instead  of  {#  , f!^}. 
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A A 

It  would  thus  appea’-  that  it  is  always  better  to  use  {'i^ , C^} 

A A 

rather  than  {.»  , c.,}*  It  should  be  remarked,  however,  that  Southwell 

A /\ 

and  Vaisey  [1946,  p.  127]  in  general  preferred  C^}  (which  they 

A,  /*V 

called  Method  A)  to  , c^}  (which  they  called  Method  B),  because 

equation  (101  provides  less  freedom  than  equation  (11)  in  determining 
i (k+11 

It  should  perhaps  also  be  remarked  that  the  choice  of  /?,  c 
may  sometimes  be  based  on  physical  analogies,  as  for  example  when 
Shaw  and  Southwell  [ 1941,  p.  7]  visualized  the  solution  of  porous  flew 
problems  for  dams  in  terms  of  gas  balloons  held  down  by  shot-bags. 

The  mathematician  may  smile  at  this  analogy  (which  was  obviously 
influenced  by  the  wartime  conditions  then  prevailing),  but  the  lack  of 
more  precise  criteria  makes  the  use  of  such  analogies  as  good  as  any 
other  method. 

In  the  next  two  subsections  two  methods  for  choosing  boundary 
conditions  are  discussed:  a method  due  to  Garabedian  [19  56]  of 

A A 

constructing  the  boundary  conditions  8,  c , so  as  to  ensure  quadratic 
convergence  with  respect  to  changes  in  the  trial  PBS;  and  a method 
due  to  Cryer  [ 1968,  197  0a]  and  McCorquodale  and  Li  [ 1971]  of  using 
integral  relations  to  derive  alternative  boundary  conditions  for  problems 
where  the  boundary  conditions  on  the  FB  involve  an  unknown  constant 
(as  is  usually  the  case  when  Bernoulli's  equation  is  used  in  fluid 
mechanics  FBPS). 


J 


In  conclusion  it  must  be  emphasized  that  although  the  above 
remarks  have  been  concerned  with  one-domain  FBPS  for  second  order 
elliptic  equations,  the  question  of  the  choice  of  boundary  conditions 
arises  in  every  trial  free  boundary  method.  We  have  not  discussed 
this  question  for  other  types  of  FBPS  because  we  are  not  aware  of  any 
work  on  the  subject. 


3.  1.  1. 


Choice  of  boundarv  condiiioris:  Garabcdian's  method 


We  now  describe  an  ingenious  way  of  choosing  B which  is  due  to 


Garabedian  [ 1956]  . Recall  that  in  a trial  free  boundary  method  the  condition 
Bx'  "■*  - 0 is  first  satisfied  on  Then  the  FB  is  moved  so  that 

the  condition  r»u  ^ = 0 is  satisfied,  but  of  course  the  condition  pu  = 0 
is  no  longer  satisfied  on  r^+^. 

A A 

The  idea  of  Garabedian  is  to  construct  B so  that  is  insensitive 
to  movements  of  r.  In  the  language  of  Garabedian,  B is  constructed  so 

A. 

that  Bu  is  stationary  with  respect  to  normal  displacements  of  the  FB.  More 

A 

precisely,  B is  constructed  so  that 


(9u)  = 0,  on  T 
n 


(1) 


where  u is  the  solution  of  the  FBP.  Garabedian's  method  may  also  be 
thought  of  as  Newton's  method  applied  to  FBPS  (Garabedian  [ 195  6a,  p.  613] 
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To  construct  « we  begin  by  observing  that  we  may  choose  the 

coefficients  e of  equation  (7)  of  section  3.1  so  that  on  r 
ij 

^ -s 

6 and  r.  are  of  the  fcrm 


&u  = (ui  - Vj_j)  + T(u  " V12)  = 


P U = U " Y1Z  = °’ 


where  r is  an  arbitrary  function.  Thus, 


= 0 r-  (un  - Yu)n  + T(U  - YI?)n  + 


12  n n 


'12 


= u - (y,,)  + t\  u - (y.  ) ] 

nn  1 11  n n 12  nJ 


since  u - y 2 = 0 on  F.  Condition  (1)  will  thus  be  satisfied  if 


-u  + (y.,) 
nn  11  n 

u - (y.J 
n 12  n 


~Unn  + (Vn 
Y11  ’ (Yl2’n 


(2) 


Let  _t  be  the  unit  tangent  to  r,  let  s denote  distance  along  F, 
and  let  k denote  the  curvature  of  r.  Then  on  r we  have 


Un  = YU’ 


u = u = (y, ) • 

t s '12  s 


(3) 


As  special  cases  of  the  Frenet  formulas  , 


- ku  = u = (y 


-KU  - U — I Y,  ~ / 5 

tt  n ss  T 12  ss 


u + ku  = u = (y,,) 
nt  t ns  "Us 


(4) 
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Finally,  the  equation  ^ \x  - 0 can  be  rewritten  in  the  form 


a'  u + 2a'  u + a'  u + b!u  + b' u = 0,  (5) 

11  nn  12  nt  22  tt  In  2 t ’ ' 


where  the  coefficients  a1  and  b1  are  linear  combinations  of  the  a 

ij  l ij 

and  b..  From  (2),  (3),  (4),  and  (5)  we  find  that 
t can  be  expressed  in  terms  of  Vjj,  (Y^,  (v^)  , (v1?) s , 

a..,  b.,  and  k.  As  an  important  special  case,  if  Yjj  and  are 

constants  and  if  3 is  the  Laplacian  then  equation  (f>)  Lakes  the  form 


u + u =0, 
nn  tt  ’ 


and  we  find  that  t = k. 

In  practice,  the  value  of  r is  not  known  because  the  values  of  k, 

Y , etc.  depend  upon  r which  is  unknown.  The  simplest  strategy  (used 

(k) 

by  Cryer  [1968,  1970a]  is  to  replace  r by  t which  is  computed  using 

(k)  (k) 

r ' instead  of  r.  Then  u is  defined  to  be  the  solution  of  the  problem. 


= 0 in  fi(k), 


S<k,u<k)  = (u<k>  - V ) + Tlk,(u(k)  - v.,)  = o oil  8D(k). 
nil  i ^ 


Cryer  [ 1968,  1970a]  found  that  the  use  of  {5’,  C } a la  Garabedian 
instead  of  {/S’,  C } improved  convergence  but  that  it  was  not  clear  whether 


convergence  was  quadratic. 
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■ . 1 . 1 . Choice  of  boundary  conditions:  integral  relations 


In  a number  of  FBPS  the  boundary  conditions  on  the  FB  involve  an 

unknown  constant  X which  must  be  determined  as  part  of  the  solution. 

For  example,  in  fluid  mechanics  FBPS  one  often  has  the  boundary  condition: 

C : u \ on  r . 
n 

( 1C  ) 

One  approach  which  has  often  been  used  is  to  compute  X by  averaging: 

. (k)  (k) 

X ~ average  u^ 


We  believe  that  this  is  not  a good  approach,  and  in  this  section  we 
indicate  how  the  method  of  integral  relations  can  be  used  to  obtain 
useful  expressions  for  X. 

We  illustrate  the  idea  by  considering  the  problem  of  a plane  jet  of  an 
ideal  weightless  fluid  flowing  from  an  orifice  (Figure  1),  following  Cryer 
[ 1968,  p.  64]-  terms  of  the  stream  function  o the  governing  equation 


and  the  boundary  conditions  are: 


M 

O 

on 

AB 

4>  - AF, 

on 

CDET 

= y, 

on 

AF, 

v.  * = Ky> 

on 

BC, 

Igrnd  v 1 - k , 

on 

CD  , 

where  V is  an  unknown  constant. 


vo  = 1 


1.  Plane  jet  flowing  from  an  orifice 


Poc 


The  pressure  and  velocity  of  the  fluid  will  be  denoted  by  p and  v 
respectively.  The  limits  of  the  fluid  pressure  and  fluid  velocity  far 
upstream  and  far  downstream  are  denoted  by  p^,  p^  , v , and  v 
respectively. 

Since  (Milne-Thompson  [ 1968 , p.  114]  ), 

v = | grad  o>  I , 

it  follows  that  far  upstream  the  fluid  flows  with  velocity  vfl  = 1 parallel 
to  the  x-axis,  while  far  downstream  the  fluid  flows  with  velocity  v = \ 
parallel  to  the  x-axis. 

The  method  of  integral  relations  involves  multiplying  the  governing 
differential  equation  by  one  or  more  functions  and  integrating  over  the 

domain  £ . The  integral  is  then  manipulated  using  integration  by  parts 


together  with  the  boundary  conditions.  In  this  way  a number  of  identities 
satisfied  by  the  solution  can  be  obtained.  These  identities  often  have 
a physical  meaning  such  as  conservation  of  mass  or  momentum,  and  this 
is  the  case  in  the  present  problem. 

From  the  principle  of  conservation  of  mass  (Milne-Thompson  [ 1968, 
P.  72]  I, 


AF  • v„  - BC  • v 
0 so 


Then,  since  CDEF  is  a streamline,  it  follows  from  Bernoulli’s  equation 
(Milne -Thompson  [ 1968  , p.  10]  ) that 


12  1 2 „„„ 
-pv  + P " 2 p V*>  P3r  °n  CDEF  ’ 


where  p is  the  density  of  the  fluid. 
Hence , 


1 2 

2 pV0  ' P0 


1 2 
2pVo= 


From  (li  and  the  boundary  condition  (!  it  follows  that  v v^ 


on 


CD.  Thus,  from  (3), 


p = p , on  CD 
00  ’ 


Next,  we  apply  Euler's  momentum  theorem  (Milne-Thompson  [1968, 
p.  79  Jl  to  the  curve  ABCDEFA: 

Resultant  thrust  on  ABCDEFA  in  the  positive  x-direction 


(2) 


(3) 


(4 1 


(51 


1 
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D E 

AF  • PQ  - J pn^ds  - J p ds  - BC  • 

rate  of  flow  of  momentum  in  the  x-direction  ) (6) 

outwards  across  ABCDEFA 

2 2 
pv^,BC  - pvQAF  , 

where  n^  is  the  component  in  the  x-direction  of  the  unit  outward  normal 
on  CD. 

Noting  (3)  and  4)  and  remembering  that  n(ds  = dy,  we  see  that,  equation  b) 
is  equivalent  to 


(AF  - DE)(pQ  - p^)  - * P vqDE  + ~ p J v2ds  pv^BC  - pv^AF  . 


(7) 


Finally,  from  (7),  (2',  and  (4),  it  follows  that 

2 E 2 2 
(AF  - DE)v^_  - 2vqv^AF  + J v ds  + vQAF  0 . 

Solving  the  quadratic  equation  (8)  for  v we  find  that 


(8) 


v AF  + [AF  • DEv  - (AF  - DE)  J v ds] 


2,^  1/2 


D 


AF  - DE 


(91 


Then,  from  (21, 


AF  : BC  v«/vo  • 


(10) 


In  solving  (8)  the  larger  root  was  chosen  because  the  smaller  root 
leads  to  a physically  unacceptable  solution.  For,  let  v be  the  smaller 
root.  Then,  from  (9), 
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V x < V Aiv  (AT  - DE) 

and  so , from  I 2 , 

BC  > (AF  - DEi  , 

which  is  physically  in. possible  (see  Figure  1). 

The  arguments  leading  to  equations  (9i  and  (10)  were  originally 
used  by  Borda  [ 1766]  ; a careful  presentation  of  Borda's  arguments  is 
given  by  Gilbarg  [ I960,  p.  340]  . Tor  the  case  considered  by  Borda  v 
could  be  determined  exactly. 

Rewriting  (9  ')  in  terms  of  X and  and  remembering  that  vQ  1 
in  our  case  we  find  that 

F 

L 2 1/2 

(AF  - DE)X  AF  ->-  [ AF  • DE  - (AF  - DE)  J (,  ) ds]  . (11) 

D X 

( k. ) 

Given  an  approximate  solution  we  can  use  (11)  to  :etermine  a 

(ki  (k) 

corresponding  approximation  X . We  may  expect  X to  be  a good 
approximation  since  the  right  hand  side  of  (11)  only  involves  the  values  of 

...  on  the  segment  DE  of  the  known  boundary,  and  it  is  reasonable  to  expect 

(k)  (k)  . ,_(k) 

- Uj  ' to  be  small  on  DE  even  if  iu  - is  large  on  I 

x x 

McCorquodale  and  Li  [1971]  use  the  same  idea  for  the  problem  of 
flow  through  a sluice  gate.  They  are  not  too  successful,  but  this  may 
be  because  they  adjust  using  a global  method  (see  section  3.2.  3). 
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3.2.  Movement  strategy 

The  most  difficult  aspect  of  trial  free  boundary  methods  involves 

the  movement  of  1 to  1'^'  Unfortunately,  many  authors  simply 

say  that  "the  trial  FB  was  adjusted  until  both  boundary  conditions  were 

satisfied"  and  give  no  further  details. 

Handworkers  such  as  Southwell  and  Vaisey  [ 1946]  did  not  use 

specific  rules  to  move  the  boundary.  It  would  of  course  be  possible 

to  emulate  this  by  using  a computer  interactively:  the  computer  would 

compute  tire  approximation  u|  and  the  error  Cu^  , while  the  human 

operator  would  adjust  the  boundary  manually  using  a light  pen  at  a 

console.  It  is  also  conceivable  that  the  boundary  could  be  moved  by  an 

"adaptive -learning " program  which  learnt  from  its  mistakes.  To  our 

knowledge  however,  in  computer  implementations  the  boundary  has 

always  been  moved  using  a definite  algorithm. 

In  computer  implementations  of  trial  free  boundary  methods  it 

(k) 

is  usually  convenient  to  regard  the  boundaries  r as  being  defined 

(k)  (k) 

by  a number  of  parameters,  a.  , . . . , a , say.  Possible  parameteriza - 

■ P 

lions  include: 

(i)  is  polygonal  with  p/2  vertices.  The  parameters  a[^,  . . . , 

a^‘ 1 then  represent  the  x and  y coordinates  of  the  vertices 
P 

, r<k> 

of  I 

(k) 

(i i i I is  the  curve 
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y(x) 


V .(k) 


u 

j ) 


J J 


where  the  functions  are  known. 

J 

„(k) 

(m)  i is  the  curve 


P(  (k) 

X(t)  = l a ' q (t)  , 

j 1 1 J 

y ft)  = alk  'tS.(t)  , 

p/1+1  J J 


where  the  functions  g are  known. 

j 

(k) 

(iv)  r is  defined  implicitly  by 


(k) 

f(x,  y,  a1  , 


a‘k’)  = 0 
’ P 


where  f is  a given  function. 


We  will  denote  the  curve  corresponding  to  a 


(k) 

1 : 


I (a 


(k)  (k)  (k)  (k)  , (k) 

1 , a ) or  F(a  ) where  a = (aj  , 


(k)  K 
a by 

P 

(k). 


Methods  for  moving  the  boundary  fall  into  three  categories  which 
we  call  local,  integral,  and  global , respectively.  To  explain  these 
methods  we  consider  the  problem  of  porous  flow  through  a dam  (see  section  0). 
In  terms  of  the  velocity  potential  u = $ the  FBP  is: 

C7u  = u + u = 0 , in  & , 

#u  = u + Ky  = 0,  on  r , 

(!u  = u - 0 , on  r t 

together  with  appropriate  boundary  conditions  on  the  fixed  boundary. 
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f\ 


Figure  1.  Mnvina  the  boundary 


rf(k  ♦ / \ 


r 


D 


/t) 


C 


(k) 

In  a local  method  the  adjustments  to  r are  made  at  individual 
points  on  the  basis  of  the  error  Cu^1  at  these  points.  Thus  Cu^',  or 

rather  Cu^\  is  computed  at  m points  ' e r ^ ; l<j<m.  If 

(k)  (k)  (k+1) 

Cu^  (P.  ‘ ) * 0 then  a nearby  point  P.  is  determined  so  that 

CuP^(p!"  **)  = 0;  that  is  p|k*  is  "moved"  to  p|k+^  (see  Figure  1).  The 
b j ’ ] j 

curve  F k +1^  = r(a(k+1))  is  obtained  by  fitting  a curve  of  the  form  T(a) 

(k+1)  k+1 

through  the  points  P.  ' . If  m - p then  l'(a^  ) will  in  general 

pass  through  all  the  points  P^k+  \ If  m > p then  F(a*+^)  is  found 

J 

by  an  appropriate  curve-fitting  procedure  such  as  least  squares. 
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is  found  by  "integrating 


In  an  integral  method  the  curve  1 


. (k  4 1 ) 


the  boundary  condition  eu  In  tho  present  example,  u 


(k 


is  a 


velocity  potential  from  which  the  velocity  field  v can  be  found  by 

differentiation.  The  condition  (>u  = u^  = 0 corresponds  to  the  physical 

condition  that  r is  a streamline.  Starting  from  the  fixed  point  A, 

the  streamline  passing  through  A can  be  determined  by  integrating  the 

velocity  field  v , cind  this  streamline  is  teken  to  dg  • 

In  a global  method  a set  of  p perturbed  boundaries 

r(kJ)  = (a[k) ajkJ,  ajk)  4 6a  j10,  a‘kj,  ...,  ajj0)  and  the 

corresponding  solutions  uk  are  generated.  This  information  makes 

(k) 

it  possible  to  estimaiL.  the  dependence  of  cu  upon  the  parameters 
a^k).  The  new  approximate  FB  is  chosen  so  as  to  minimize 

the  error  cV^1*.  In  particular,  if  it  is  assumed  that  Cu(k'  depends 
linearly  upon  a(k),  then  a(k  1]  can  in  general  be  chosen  so  that 

(k+1)  . . . . (k+1) 

Cu  is  zero  at  p points  on  i 


These  three  methods  of  moving  T are  discussed  in  greater  detail 
in  the  following  subsections.  In  the  final  subsection  we  discuss  the  case 
when  the  boundary  conditions  involve  an  unknown  constant. 
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3.2.1.  Movement  strategy:  local 

(k) 

In  a local  method  of  moving  the  boundary  the  error  cvi^  in  computed 

at  m points  p(k)(](k),  l<j<m.  11  cu!k)(p|k))  i 0 then  a nearby  point 
j n j 

P<k is  determined  so  that  the  boundary  conditions  are  "satisfied  better  ' at 
j 

P(k,])-  that  is  P(k)  is  "moved"  to  P(k+1>.  The  curve  r(k_tl)  is  found  by 
j ’ ’ j J 

fitting  a curve  through  the  points  P 

Handworkers  such  as  Southwell  and  Vaiscy  [ 1946]  did  not  use 
specific  rules  to  move  the  boundary  and  used  trial  and  error.  Indeed,  when 
comparing  the  choice  of  a Dirichlet  boundary  condition  for  C (which  they 
call  Method  B) 


V=Un+  nr0’ 


ciu  = u + nz = °’ 


(i) 


with  the  choice  of  a Neuman  condition  for  C (which  they  call  method  A) 


92u  -- u + v21  = 0, 


(2) 


C u = u +\  =0, 

2 n 22 


Southwell  and  Vaisey  [ 19  46,  p.  127]  comment  that  they  prefer  Method  A to 
Method  B because 

"Method  B . . . yields  a definite  indication  regarding  the 
shape  to  be  adopted  in  the  next  stage  of  the  computations 
...  in  Method  B divergence  is  not  under  control". 
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Southwell  [1946,  p.  Z25]  summarizes  his  viewpoint  very  well: 

"The  ability  of  Relaxation  Methods  to  deal  with  problems 
such  as  these  is  due  to  the  retention,  throughout  their 
development,  of  a tentative  Quality  which  orthodox 
methods  do  not  possess.  There  is  nothing  new  to  the 
notion  of  continued  approximation  to  a wanted  solution, 
for  this  is  the  basis  of  all  'iterative'  methods  of  attack. 

What  is  novel  is  the  freedom  that  throughout  is  left  to 
the  computer,  to  decide  the  nature  of  his  next  step.  " 

It  is  of  interest  to  mention  some  of  the  details  of  trial -and -error 
methods : 

(i)  Vaisey  [ 1961]  has  informed  us  that  in  the  computations  the  ratio 
(local  correction)/ (local  error)  sometimes  changed  sign  during 
the  process  of  solution.  Furthermore,  in  some  cases  the  error 
over  a section  of  the' boundary  could  only  be  reduced  by  adjusting 
a different  section  of  the  boundary.  The  problem  of  a waterfall 
(Southwell  and  Vaisey  [1946,  p.  144])  presented  considerable 
difficulties,  and  we  have  been  told  that  in  some  cases  it  was  not 
possible  to  reduce  the  error  by  moving  the  boundary  of  the  jet 
and  that  instead  the  edge  of  the  waterfall  was  moved. 

(ii)  Abul-Tetouh  [ 1949]  and  B.  W.  Hunt  [ 1967]  considered  the  problem  of 
an  axially  symmetric  jet  from  an  orifice.  Following  Southwell 

and  Vaisey  [ 1946,  p.  134]  they  used  the  Neumann  boundary  condition 
for  c which  takes  the  form 


9 


1 


- \ 


x an 


(3) 


-85- 


r 


Abul - I'otouh  u.-.i-d  finite  dilfcrences  and  IIu.it  used  an  integral 
equation  to  compute  uj  \ but  botii  authors  used  trial -and-error 
to  move  the  I B and  give  clear  descriptions  of  their  experience. 


Abul-Fotouh  [ 19 49  , P-  25]  writes: 

"The  v due  of  g indicates  the  difference  between  the  asymptotic 
velocity  and  the  local  velocity  at  a certain  point  on  the  surface, 
and  this  should  be  reduced  to  zero.  The  tendency  of  the  variation 
of  g will  determine  the  new  assumed  changes  in  the  curve.  A 
positive  value  of  g shows  that  the  local  velocity  is  too  large, 
and  the  jet  section  should  increase.  The  curvature  of  the  stream 
surface  and  the  cross  section  of  the  flow  affect  the  magnitude  of 
the  velocity  at  any  point  on  that  surface.  An  increase  in  cross 
section  decreases  g and  a decrease  in  curvature  decreases  g 
also.  These  two  factors  will  work  together  to  produce  the  final  g. 
The  method  of  choosing  the  right  curve  to  make  g approach  zero 
is  not  easy;  it  requires  considerable  experience  and  personal 
judgment.  However,  some  guidance  is  at  hand  if  one  knows  that 
an  increase  in  the  jet  radius  at  any  point  will  decrease  the  value 
of  g at  that  point  and  increase  the  values  of  g at  the  adjacent 
points  on  both  sides  of  the  considered  point,  since  the  curve  be- 
comes more  concave.  If  the  values  of  g alternate  in  sign,  the 
shape  of  the  curve  should  change.  If  g has  the  same  sign,  the 
asymptotic  radius  should  change  with  the  corresponding  change 
in  the  curve.  Trial  and  error  will  show  the  tendency  and  the 
amount  of  change  in  the  ordinates  of  the  curve.  No  definite  value 
of  the  change  of  the  ordinates  can  be  expressed  to  produce  a 
certain  change  in  the  surface  velocity  at  a certain  point,  since 
the  curvature  of  the  profile  affects  the  results  to  some  extent. 
However,  as  a rough  estimate,  a change  of  ^ 0.004  in  the  ordinate 
at  a point  will  change  the  velocity  at  that  boundary  point  by  ± 1 ro. 
Near  the  orifice  edge,  the  situation  is  more  difficult,  since  the 
profile  curves  sharply  in  this  zone  and  it  is  only  by  trial  and  error 
that  one  can  arrive  at  the  right  shape.  If  two  assumed  curves 
happen  to  produce  different  signs  of  g for  similar  points,  then 
the  right  curve  lies  between  these  two  and  a simple  interpolation 
for  the  ordinate  will  be  of  great  help. 


J 
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Hunt  | 1967 , p.  18]  writes: 

"Details  of  the  method  used  for  locating  the  free  stieamline 
consisted  of  initially  choosing  an  asymptotic  radius  r^.  Next, 
the  shape  of  the  free  streamline  between  the  orifice  and  the 
region  of  uniform  flow  was  adjusted  by  trial  and  error  until 
velocities  along  the  portion  of  the  free  surface  farthest  from 
tlic  orifice  were  essentially  constant,  linally,  if  free-surface 
velocities  in  the  vicinity  of  tire  orifice  were  higher  than  the 
chosen  asymptotic  velocity,  the  asymptotic  radius  was  decreased. 
Conversely,  if  velocities  in  this  region  were  too  low,  the  asymp- 
totic radius  was  increased.  The  entire  process  had  to  be  repeated 
a number  of  times  until  a satisfactory  free- streamline  geometry 
was  found.  This  usually  meant  calculating  velocities  along  the 
free  surface  for  ten  to  twenty  different  geometries  in  order  to 
obtain  one  solution. 


In  addition  to  the  authors  mentioned  already,  the  following  authors 
adjusted  the  FB  by  trial-and-error: 

With  finite  differences:  Vandrey  [ 1940]  , Shaw  and  Southwell  [ 19  41]  , 

Binnie  and  Davidson  [ 1949]  , Yang  [ 1949]  , Rouse  and  Abul-Fetouh 
[ 19  50]  , van  Deemter  [ 19  50]  , Ball  [ 19  51]  , Boulton  [ 19  51]  , 

Brunauer  [ 19  51]  , Citrini  [ 19  51]  , Kashef,  Toulouklan  and  Fadum 
[ 19  52]  , Allen  and  Dennis  [ 19  53]  , Boreli  [ 19  55]  , H . P . Hall  [ 19  53]  , 
McNown,  Hsu,  and  Yih  [ 1955]  , Young,  Gates,  Arms,  and  Eliezer 
[ 1955],  Dumitrescu,  Ionescu,  and  Toth  [ 19  56]  , Jakobsson  and 
Floberg  [ 19  57]  , J.  A.  Murray  [ I960]  . 

With  integral  equations:  Trefftz  [ 1914,  1916]  , Wagner  [ 1932]  , 


Schach  [ 19  3 5]  . 


With  the  advent  of  computers  it  became  desirable  to  automat 


the  movement  of  the  FB.  The  approach  usually  followed  has 

D(ktl) 


determine  the  points  P 


1 


according  to  the  condition  that 


C!u(k)(P(h+1))  = 0 

n ] 


together  with  the  condition  that 


(k+ 1 ) (k)  (k)  (k) 

P.  - P.  = a d 
) ) ) ~j 


where  d|k*  is  a specified  unit  vector. 

(k) 

Possible  choices  for  d,  are: 

~) 

(k)  (k) 

(a)  The  unit  outward  norrrial  to  r at  P. 

(b)  The  unit  vector  in  one  of  the  coordinate  directions. 

(k)  (k) 

(c)  The  unit  outward  conormal  to  r at  P . That  is 


governing  equation  is 


C7u  = (f.u  ) + (f  u ) 

1 x x 2 y y 


and  the  unit  outward  normal  is 


n = (n^  n2), 


then  the  unit  outward  conormal  is  the  vector 


- = (fini>  f2n2} A (flnl)2  + (f2n2)2^/2 


been  to 


(4) 


(5) 


, if  the 
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(The  advantage  of  u.  mg  the  conormal  is  that  t) i o boundary  conditions 
for  £7  arc  often  expressed  in  terms  of  the  conormal.). 

(d)  To  specify  d.  as  part  of  the  input  data  as  is  done  by  R.  L.  Taylor 
[ 1966]  (see  also  Kenly  and  Busch  [ 1971]  ). 


The  computation  of 
upon  the  structure  of  r . 


(k4  1) 

P;  so  as  to  satisfy  (4)  depends. 

) 

The  most  obvious  approach  (to  us)  is: 


of 


course, 


Method  1:  Define 


compute  or  estimate 


,(k)  (k)  D(k)  x ,(k) 

f (a)  = c u (P.  + acb  ); 


f f,k,M 

d«  ; 


tt=0 


and  then  set 


(k)  f(k>,  \ I , 

a.  - -f.  (0  /(—  f.  (a)  n) 

j j aa  j \a- 0 


(6) 


However,  various  other  approaches  have  been  used: 


Method  2:  If  Cu(P)  is  of  the  form 


Cu(P)  ■ F (u(P) , u (P),  P)  = 0 


then  some  authors  define  a\^  by  means  of 

) 
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(7) 


r(uV!'V  u(kV::,i.  P<k,+  o!,t,d!k,,  = 0: 

j n j j J -j 


that  is  the  dependence  of  u*k>(Pi  upon  P is  neglected. 


Method  An  even  cruder  approach  is  to  introduce  a function  G(x) 

(such  is  G(x)  x)  and  set 


cokl  G(ru(P(k)))  . (8) 

J J 

Although  Method  1 is  more  elegant,  Methods  2 and  3 have  been  used  with 
great  success  in  various  problems,  and  it  is  conceivable  that  in  some 
cases  Methods  2 and  3 are  more  stable  than  Method  1. 

We  now  illustrate  the  above  remarks  by  some  specific  examples. 
Most  authors  have  chosen  C to  be  the  Dirichlet  condition 


Cu(P)  - u(P)  + \12(P)  = 0 . (9) 

With  this  choice  of  e the  approaches  which  have  been  used  include. 

(i)  In  the  programs  of  R.  L.  Taylor  for  porous  flow  problems  which  have 
been  used  by  Taylor  and  Brown  [ 1967]  , Kealy  and  Busch  [ 19  < 1]  , 
Kealy  and  Williams  [ 1971]  , and  Pettibone  and  Kealy  [ 1971]  , the 
boundary  condition  is  of  the  form 

C : u = 0 

where  u is  the  pressure.  Taylor  (see  R.  L.  Taylor  [ 1966 , p.  6] 

and  Kealy  and  Williams  [1971,  p.  62])  determines  P^'11  by 

p(Hl)  (kL  pJk)  p(kl|tl  (ki 

j j h J -J 
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where 

p is  a cnrrcc 

tion 

tactor.  (This  is  M ,-thod  i with  G(x)  Px). 

both 

i . i nd  t he  dii u 

•••tii 

d^'1  must  be  provided  by  file  user  as 
-j 

data. 

R.  L.  Tuylor  [ 

19  06 

] suggests  the  choice  (1  . ZS  while 

Ke-ily 

and  Williams  | 

1 1971 

) suggest  p - .7.  No  advice  is  given  about 

(k ) 

the  choice  of  the  directions  d.  but  they  should  obviously 
make  an  angle  of  less  than  90°  with  the  positive  y-axis. 

This  is  a very  crude  approach  - but  it  works. 

In  porous  flow  problems  the  FB  condition  can  also  be  written 
in  the  form  (see  section  0.) 


C : u - y = 0 


. „ D(k+1) 

where  u is  the  hydraulic  head.  This  suggests  that  P 
(X|k+1)(  y|k+1))  be  defined  by  moving  P(k)  in  the  direction  of 

the  y-axis , 


(k+1)  (k) 

x.  = x.  , 
J ) 


(k+1) 


yi 


u^(P3(k>) 


(10) 


(Method  2)  and  this  has  been  used  by  Neuman  and  Witherspoon 
[ 1970]  . Neuman  and  Witherspoon  also  incorporate  two  other 
features:  a preliminary  calculation  is  made  to  determine  the  flow 
on  the  seepage  surface;  and  certain  restrictions  (which  they  describe 

D(k+ 1) 

in  some  detail)  arc  imposed  on  1. 

Fenton  f 197  2a,  p.  4-4  and  Appendix  D]  (see  also  Parkin  j lc>7l]  ) 
replaces  (H)  with 
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(11) 


(k+1)  (k)  . | (k)  (k)  (k) 

x,  = 4 d^luh  (P|  ) - y,  ], 


k+1  k , , (k).p(k)  (k) 

y = y + «LU.  (P.  - y . , 

y x h j j 


where  is  taken  to  bo  1.  5 and  where  the  constants  are  chosen  so  that 

P(k1’  lies  on  a certain  line.  (This  is  Method  3 with  G(x)  - ^x). 

j 

(iii ) [’or  the  Dirichlct  condition  (9)  Method  1 may  be  written  in  the  form 


(k) 


4 Vl7,P,k>, 

J U ..J 


7^  [u(k,(P)  + y1?(P)]  (k) 

12  p-  pv  ' 


(12) 


J 

The  quantity  9^  /a»  is  readily  computed.  If  the  first 

boundary  condition  is  chosen  to  be  of  the  form 


#u=|^+y  =0, 

9o  22 


then 


9u 


(k) 


9 a 


“V 


22 


(k) 


otherwise  then  — can  be  estimated  numerically  or  otherwise. 

9 £ 

This  approach  was  used  by  Garabedian  [ 19 56 , 1956a],  Cryer[  1962, 

1968,  1970a]  and  J.  M.  Taylor  [1971,  p.  35]  (see  also  Aitchison  [ 197  2] ). 

The  same  approach  was  used  in  a special  case  by  Tu  [1971,  p.  43]  who 

(k) 

introduced  a relaxation  parameter  w so  that  a was  determined  by 

u(k)(P(k>)  + y,  (P|k)) 

3 * ^ J 


(k)  = . 

' “ ^ [u(k)(P)4V,  ,(P>] 


9« 


12^'VpOo 
" i 
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In  a number  of  cases  the  boundary  condition  C ha."  boon  taken 
to  be  in  the  Neumann  form 


eu(P)  un(P)  •<  y (P)  o. 


With  this  choice  of  C the  approaches  which  have  been  used  include: 


(i)  In  fluid  mechanics  TBPS  involving  gravity  we  have 


Cu  = (~)2  - 2gy  = 0, 
a n 


so  that  we  may  set  (Method  2) 


(k+1)  (k) 

x.  = x. 

J J 


(k+l)'_  f|u  (p(k)  Z,  _ 

) an  j 


Stewart  [ 196  5,  p.  42]  and  Stewart  and  Davidson  [ 196?]  have  used  this  method  , 
Finnemore  and  Perry  [ 1968]  consider  a porous  flow  problem  and 
write  the  boundary  condition  in  the  form 


au  . zn 

Cu  = — - sin  0 

ay 


where  0 is  the  angle  between  the  normal  to  r and  the  x-axis. 

(k)  (k+1)  (k) 

Finnemore  and  Perry  move  P]  by  setting  x]  = x.  ' and 

setting  y!^+^  equal  to  a rather  complicated  and  arbitrary  fur  w 
V fk ) 

of  Cu  (Pj  (Method  3). 
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In  a number  of  instances  it  has  been  found  desirable  to  smooth 
(k 1 1 ) 

the  points  Pj  so  as  to  prevent  undesirable  oscillations: 

(a)  Fi nnemorc  and  Perry  [ 1968]  used  a standard  smoothing  subroutine. 

(b)  Cryer[1962,  1968,  1970a]  determined  m points  p|k  \ 1 < J < m, 

and  then  fitted  a curve  r (a ) of  prescribed  form  through  these 
_(k+l)  ,,,  ktl 

points  to  obtain  l = Ma.  ). 

Stevens  [ 1974]  solves  an  axisymmetric  magnetohydrostatic  FBP 
and  apparently  uses  a local  method  (he  gives  few  details). 

In  summary,  local  methods  have  been  used  very  successfully 
despite  their  apparent  arbitrariness. 


3.  2.  2.  Movement  strategy:  integral 

In  an  integral  method  of  moving  the  I'B  the  boundary  condition 
Cu  : 0 is  expressed  in  an  implicit  form  such  as 


F (u(x) , u (s),  u (s),  x(s),  x(s))  = 0, 
’ x ’ y 


(1) 


where  the  FB  is  the  curve  x = x(s)  and  where  x(s)  denotes  the  derivatives 
of  x(s).  Given  I'(k)  and  u(k),  the  curve  r(k+1)  is  obtained  by  integrating 
approximately  the  differential  equation  for  x(s), 


F(u[]°(s),  u^s),  U-,  x(k4l)(s),  x(k+1)(s))  = o.  (2) 


(s) 


In  (2)  uf^  (s)  etc.  represent  approximations  to  the  values  of 
etc.  at  the  point  x^k  *\s).  If  x^  *(s)«j9  then  u ^ (s)  can 
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be  c iitai  ? i by  interpolation . If  x^  ^(s)/ & then  (s)  must  be 
obtained  by  extrap<  dation . In  the  applications  the  method  used  to  compute 

u^'(s)  does  not  seem  to  be  critical;  often,  u ' (s)  is  taken  to  be  the  value 

(k)  (k+1),  . 

of  u at  the  gridpoint  nearest  to  x (s). 

The  following  applications  of  integral  methods  have  been  made: 

(i)  The  self-consistent  method  of  computing  the  magnetopause  which 

was  introduced  by  Beard  in  I960  and  implemented  numerically  by 

Mead  and  Board  [ 1964]  , Baker,  Board,  and  Young  [ 1964]  , is  an 

integral  method.  In  the  general  case  considered  by  Olson  [ 1963, 

1969]  the  I'B  is  expressed  in  the  form 

r = R(0,  4>) 

where  (r,  0,  $)  arc  polar  coordinates.  The  expression  for  Cu  s F 
is  given  in  section  2.4.3. 

(ii)  von  Kerczek  [ 196  5]  used  a version  of  the  integral  method  for  two 
problems,  the  flow  over  a submerged  vortex  (successfully)  and 
flow  under  a planing  surface  (unsuccessfully).  The  method  of 
von  Kerczek  differed  from  the  integral  method  as  described  above 
in  that  during  each  iteration  only  one  boundary  point  was  moved. 

(iii)  Mogal  and  Street  [ 197  2,  197  4]  solve  the  problem  of  a plane 
Riabouchinsky  cavity.  For  this  problem, 


Cu  = F = ^ x(s)  + ib  y(s) 
x y 


where  x = (x,  y). 
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(iv)  Tanner  [ 197  3]  and  Nickell,  Tanner,  and  Caswell  [ 197  4]  compute 
the  flow  of  a viscous  jet  from  a tube. 

(v)  Comincioli,  Guerri,  and  Volpi  [1971]  and  Baiocchi,  Comincioli,  Guerri, 
and  Volpi  [ 1973]  solve  the  problem  of  porous  flow  through  a dam. 

(vi)  Michael  and  O'Neill  [ 197  2a]  solve  a problem  involving  a fluid 
film  in  an  electric  field.  This  is  an  interesting  problem  in  that 
the  equation  (2)  for  x(s)  is  a second  order  differential  equation 
which  gives  rise  to  a two-point  boundary  value  problem. 

As  the  above  examples  show,  the  integral  method  has  proved  of 
great  value,  and  it  seems  clear  that  many  more  applications  will  be  found. 


3.2.3.  Movement  strategy:  global 

iy  \ 

In  a global  method  of  moving  r"  a set  of  perturbed  boundaries 

r<M>.r(.*> .«],  .J» 4k>> 

and  corresponding  solutions  u^’^  are  generated.  This  information  makes 

(k) 

it  possible  to  estimate  the  dependence  of  the  error  eu  upon  the  parameters 

(k'll  ...  ..  (k+1) 

a . The  new  TB  1 is  chosen  so  as  to  minimize  the  error  cu 

In  order  to  minimize  the  error  qu^  ^ we  must  have  a measure  for 
this  error,  F(a)  say.  Two  choices  for  F have  been  used: 

(i)  F1(a(k))  - (Flj(a(k)))  = (Cu£k)(pJk)>>; 

that  is,  F^a^')  is  the  m-vector  of  the  errors  at  the  m points  pjk). 
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M 


cm 

j=l  ‘ 

that  is,  Tz  is  the  least  squares  error. 


(k) 


When  using  Fj  it  is  assumed  that  depends  linearly  upon  a : 


F (a,kl  ♦ 6a)  = A,k,6a  * b‘k> 


ty  \ 

The  p x p - matrix  A and  the  p-vector  b ' are  found  by  computing  p 


perturbed  solutions 


u(k’j\  j = 1,  . . . , p.  Then 


(k+1)  (k)  /A(k)  -1  (k) 

a = a - (A  ) b 


This  approach  has  been  used  by  Sankar  [1967,  p.  193]  and  Fox  and  Sankar 
[ 197  3]  for  tiio  problem  of  an  axially  symmetric  Riabouchinsky  cavity.  Tu 
[1971,  p.  60]  used  the  same  method  for  the  problem  of  flow  over  a circular- 
crested  weir,  but  Tu  also  intervened  manually. 

Midgley  [ 1963]  and  Midgley  and  Davis  [1962,  1963]  use  a similar 
method  for  computing  the  magnetosphere  in  a plasma  of  constant  pressure, 
and  in  a solar  wind,  the  difference  being  that  the  components  of  F^  are 
taken  to  be  the  moments  of  the  surface  current  (see  section  2.4.31. 

Although  the  choice  F^  (least  squares)  seems  reasonable,  it  has 
proved  to  be  unsatisfactory.  Arms  and  Gates  [ 19  57]  were  unable  to  obtain 
convergence  for  the  problem  of  a Riabouchinsky  cavity.  We  used  least 
squares  on  a problem  arising  in  the  theory  of  stellar  evolution  (Cryer  [ 196  2, 
chapter  10,  p.  7])  and  also  failed  to  obtain  convergence.  Rippin  [ 1959] 
and  Rippin  and  Davidson  [ 1967]  used  least  squares  for  the  problem  of  an 
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axially  symmetric  spherical  cap  bubble,  but  had  to  intervene  manually.  The 

(k) 

difficulties  with  the  choice  F^  appear  to  be  due  in  part  to  the  fact  that  Cu  can 
only  be  computed  to  a few  decimal  places  together  with  the  fact  that  least 
square;  problems  lire  ill-conditioned  unless  care  is  taken.  McCorquodale 
and  Li  | 19  71]  use  least  squares  for  the  problem  of  flow  under  a bevelled  sluice 
gate.  McCorquodale  and  Li  [ 19  71,  Figure  3]  plot  F^  as  a function  of  their 
parameters  a^  and  a^  and  this  illustrates  a difficulty  with  least  squares 
which  we  also  experienced:  the  minimum  of  is  not  clearly  defined  and 
It  is  not  clear  that  there  is  a point  at  which  T = 0. 

We  conclude  with  some  remarks: 

(a)  Whatever  the  choice  foi  F,  it  is  clear  that  it  is  very  desirable 

to  choose  the  perturbations  so  that  they  are  orthogonal  or  nearly  so. 

(b)  It  is  necessary  to  ensure  that  F depends  smoothly  upon  a. 

Arms  and  Gates  [ 1957,  p.  6]  used  finite  differences  with  a fixed 

(k) 

grid.  The  perturbations  of  r sometimes  resulted  in  the  FB 

moving  across  a gridpoint,  and  it  was  found  that  this  caused 

, , , fk) 

unpleasantly  large  changes  in  qu 

(c)  In  comparing  global  methods  with  local  and  integral  methods,  it 
must  be  remembered  that  each  step  of  a global  method  requires 
approximately  p times  as  much  work  as  one  step  of  a local  or 
integral  method. 

(d)  In  summary,  global  methods  of  moving  the  FB  have  not  been  very 
successful  in  comparison  with  local  and  integral  methods.  We 
suspect  that  this  is  because  global  methods  do  not  take  advantage 


of  any  properties  of  the  solution. 
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3.2.4.  Estimation  of  unknown  constant 


In  many  inv;  ;cid  Uui  1 mechanics  HIPS  one  boundary  condition  on 
the  I'B  is  the  Bernoulli  equation 

~ Iqrad  u |Z  4 pqy  - \ = 0, 

where  X.  is  a constant.  In  a few  problems,  such  as  the  problem  of  a 
Borda  mouthpiece  the  value  of  X can  be  determined 
a-priori.  In  most  problems,  however,  it  is  necessary  to  determine  k as 
part  of  the  solution.  It  is  of  interest  that  porous  flow  TBPS  with  sheet  piles 
also  have  an  unknown  constant  in  the  boundary  conditions  (see  Cryer  [1976, 
section  3.1.3]). 

Methods  which  have  been  used  to  determine  the  unknown  constant 
X.  include: 

(a)  Using  a boundary  condition  B which  does  not  involve  \ to  compute 
u^,  evaluating 

MP(k))  = - I grad  u(k)(P(k))  I2  + pgy|k>  , 
j 2 J 3 


and  setting 


. (k+1)  WP(M. 

X - average  MF.  ) . 


Using  integral  identities  as  described  in  section  3.1.2. 


A 


3.3.  A convergence  proof  for  a model  problem 

So  for  as  wc  are  aware  no  proof  of  convergence  of  a trial  free 

boundary  method  has  been  given.  Cryer  [ 1968,  1970a]  analysed  a very 
simple  model  problem  and  for  lack  of  better  results  this  work  is  described 

below . 

The  model  FBP  to  be  considered  is  as  follows:  it  is  required  to 
find  u satisfying 

u + u = 0,  in  ^ , 
xx  yy 

where  6 is  of  the  form  shown  in  Figure  1. 


The  FB  r is  the  curve  BC.  The  auxiliary  restraints  are  that  r should 
pass  through  the  fixed  point  C and  that  on  r y should  be  a monotone 
decreasing  function  of  x.  The  boundary  conditions  are: 
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u - ( 1 - y ) , on  AC, 


/9u  = 0 / un  •<  v 10  , on  BC, 


u - 1 on  AB, 

n ’ 


Cu  - 0 - u - ~ , on  BC. 


It  is  readily  verified  that  a solution  of  this  problem  is  given  by 


u : u = 1 - y - 3x, 


r : y = Z ' 3x 


We  believe  that  this  solution  is  unique  but  have  not  tried  to  prove  this. 

Since  the  FB  is  a straight  line  it  has  curvature  k = 0.  It  follows  that 
the  boundary  conditions  are  already  in  the  form  suggested  by  Garabedian  [ 19  5o] 
(see  section  3.1.1). 

To  solve  the  problem  the  approximate  TBS  are  assumed  to  be  straight 

(k) 

lines  passing  through  C.  That  is,  T is  assumed  to  be  of  the  form 


1 x (k) 
y ” 2 + c >< 


The  problem 


(k)  , (k)  n fi(k) 

u I u = 0 , i n & 
xx  yy 


*u(k)  - 0,  on  D*(k) 


J 


I 


ho:;  the  exact  soiutioi 


u(k)  1 y + x{|  10([c(k)]Z  4 1)]1/Z  - 1}/C(k) 


The  condition  that  Cu'k*  = 0 on  F^k+^  is 
is  defined  by 


is  satisfied  exactly  if  c 


(kil) 


c(k4l)  = {[  I0([c(k)]  2 + 1)]1/2  - l}/c(k) 


Thus  in  this  very  simple  problem  the  approximate  solutions 


(k) 


u and  th e 

approximate  FBS  1'  ’ are  known  exactly. 

To  analyse  the  behaviour  of  the  coefficients  c^  it  is  helpful 
to  observe  that  if 


f (c)  = [1  + c2]  1/2  - 101/2  , 


then 


c(k4l)  = c,k)  - f(c(kVf(dk)] 


(k) 

so  that  the  sequence  c is  identical  with  the  sequence  which  would  be 
obtained  by  starting  with  the  initial  guess  c^°*  and  applying  Newton's 
method  to  the  equation  f(c)  = 0.  Noting  that  f(c)  is  convex  for  c < 0 
and  that  f(o)  < 0 we  have  (Henrici  [ 1964,  p.  79] ): 

Theorem  1 


(k) 


Tor  any  initial  guess  c(o)  < 0 the  sequence  of  approximate  TBS 
converges  quadratically  to  the  TB  1\  FJ 
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L. 


k 


We  coiiclii'li  the  discussion  of  the  model  problem  witii  some 
remarks : 

(i)  The  fact  that  1 ^ converges  quadratically  to  I supports  the 

assertion  of  Carabedian  | 19  56,  p.  223]  that  quadratic  convergence 
is  obtained  by  choosing  & and  C according  to  his  method.  The 
fact  that  sequence  c ^ can  be  obtained  by  applying  Newton  s 
method  to  the  equation  f(c)  0 reminds  one  that  Carabedian' s 
method  can  be  thought  of  as  an  attempt  to  apply  Newton's  method 
to  FBPS  (Carabedian  [1956a,  p.  613]). 

(iii  Since  u1'"'  is  linear  it  can  be  computed  by  many  methods . Cryer 
[1968,  p.  46]  computed  u^'  using  finite  differences  and  found 
that  round-off  errors  did  not  affect  the  quadratic  convergence. 

(iii)  The  model  problem  was  chosen  because  it  resembles  a much  mere 
complicated  TBP  arising  in  the  study  of  stellar  evolution.  The 
numerical  results  for  the  more  complicated  problem  behaved 
similarly  to  those  for  the  model  problem  (Cryer  [ 1968]  ). 


Summary  of  numerical  experience 


In  this  section  we  summarize  the  available  evidence  about  trial 
free  boundary  methods.  Most  authors  are  rather  reticent  about  describing 
unsuccessful  methods  or  methods  which  required  a great  deal  of  adjustment 
before  they  worked.  This  is  understandable  but  it  is  unfortunate  for  two 
reasons:  (i)  later  workers  ore  not  forewarned  and  repeat  the  mistakes; 
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(i)i  the  devclop.nc-.it  of  a 


satisfactory  theory  is  hindered  since  such  a theory 
must  be  formulated  so  a::  to  exclude  non -convergent  examples.  V/e,  therefore, 
emphasize  reported  cases  of  non -convergence. 

When  the  evidence  is  examined  it  becomes  clear  that  the  behavior 
of  trial  free  boundary  methods  is  heavily  dependent  upon  the  type  of  problem:: 

Porous  flow  F3PS 

There  is  overwhelming  evidence  that  trial  free  boundary  methods 
are  stable  and  converge  satisfactorily  for  porous  flow  problems.  The  trial 
FBS  can  be  moved  using  any  local  method  and  smoothing  is  not  required. 

Most  authors  take  the  auxiliary  condition  to  be  the  Dirichlet  condition 
h y but  this  is  not  necessary.  We  summarize  the  evidence: 

(i)  Finite  difference  methods  were  successfully  used  by  a great  many 
hand  workers.  Gillian  Vaisey  [1961]  has  said  that  she  regarded 
the  problem  of  porous  flow  through  darns  as  an  easy  problem  and 

assigned  it  as  an  exorcise  to  her  students  for  hand  computation. 

(ii)  Wyckoff  and  Heed  [ 193  9]  and  later  workers  such  as  Childs  used 

conducting  paper  and  adjusted  the  trial  I'B  by  cutting  the  paper. 

Clearly  this  approach  requires  that  the  method  be  extremely  stable. 

(iii)  Recently  both  finite  difference  and  finite  element  methods  have  been 
used  with  automatic  adjustment  of  the  FB. 

(iv)  Porous  flow  FBPS  have  been  solved  automatically  using  resistance 
networks  by  Karplus  [ 19  96]  and  Flerbert  and  Rushton  [ 1966]  . 
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To  out  knowledge  only  one  case  of  divergence  for  a porous  flow  FBF 


has  been  reported  in  the  literature.  R.  L.  Taylor  and  Brown  j 1967,  p.  30 
in  : p.  32]  : < nti  n in  “ambiguity  effect'1  at  the  point  where  the  J'B  Intersects 
the  seepage  surface.  Taylor  and  Brown  do  not  give  details  but  Neuman  u 
Witherspoon  [ 197  0,  p.  89  5]  used  the  program  of  Taylor  for  a dam  with  a 
toe  drain  (Figure  1). 


Figure  I : Dam  with  a toe  drain 

The  successive  approximations  to  the  TB  converged  except  in  the  neighbour- 
hood of  the  tec  dam  where  the  FB  oscillated.  Neuman  and  Witherspoon  [ 1970] 
implemented  a modification  of  Taylor's  method  which  did  converge.  The 
reason  for  the  non-convergence  of  Taylor's  results  is  not  clear.  Neuman  and 
Witherspoon  [1970,  p.  891]  state  that  the  reason  for  the  non-convergence  is 
that  the  length  of  the  seepage  surface  is  not  known.  It  seems  to  us  that  the 
non-convergence  may  be  connected  with  the  peculiar  boundary  conditions 
on  an  overhanging  seepage  surface  (see  Cryer  [ 1976 , section  2.2.]). 

Kealy  and  Busch  [1971],  who  used  Taylor's  program,  also  mention 
the  ambiguity  effect. 


i 

I 

I 
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Most  authors  use  a reasonably  good  initial  guess  for  I 


I . M.  Taylor  [ 1971,  p.  37]  mentions  that  difficulties  occurred  when 
was  too  far  from  I . 


We  make  the  following  conjecture: 


Conjecture  1 

In  a porous  flow  FBP,  if  the  initial  guess  1 for  the  FB  is 

(k) 

horizontal  then  the  successive  trial  FBS  F converge  monotonely 
downwards.  (A  similar  assertion  was  made  by  F.  S.  Shaw  and  Southwell  [1941, 
p.  8]  who  asserted  that  it  could  be  proved  by  the  maximum  principle.  ) ■ m 


Fluid  mechanics  FBPS 

Trial  free  boundary  methods  for  fluid  mechanics  problems  often 

exhibit  instability.  Various  smoothing  techniques  seem  to  be  necessary: 

(k+1)  (k ) 

obtaining  the  new  approximate  FB  F by  moving  the  points  on  F 

and  then  fitting  a curve  through  these  points;  using  a weighted  average  of 

the  new  and  old  solutions;  manual  intervention. 

Southwell  and  Vaisey  [1946,  p.  127]  who  solved  eleven  fluid 

mechanics  FBPS  by  hand  observe  that  "usually  the  shape  [of  successive 

I ^1]  is  found  to  diverge  in  some  part  at  least". 

Many  authors  have  reported  convergence  difficulties  when  adjusting 

the  FB  automatically:  Young,  Gates,  Arms,  and  Eliezer  [ 1955,  p.  15]  and 

Arms  and  Gates  [ 19  57,  p.  5]  could  not  solve  the  Riabouchinsky  cavity 
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k 


problem;  von  Kerzcek  | 196s,  p.  56]  could  not  solve  the  problem  of  a flat 

plate  hydroplane:  Stewart  | 196  S,  p.  92]  was  forced  to  intervene  manually 

when  solving  a bubble  problem:  Tu  | 1971,  p.  87]  obtained  divergence  when 

solving  the  problem  of  a curved  weir  with  gravity;  Steffen  [ 1973]  found  that 

the  iterates  did  not  converge  when  solving  the  problem  of  an  axially 

symmetric  wall  jet.  We  had  personal  experience  of  such  instability  when, 

around  i960,  we  tried  to  use  a computer  to  solve  the  progressive-wave 

problem  which  had  previously  been  solved  by  Southwell  and  Vaisey  [1946]  : 

(k  i 

many  choices  of  9 and  C were  tried,  but  in  each  case  the  curves  r 
rapidly  developed  undesirable  wiggles. 

Other  TBPS 

There  is  not  sufficient  evidence  about  other  TBPS  to  draw  any  definite 
conclusions,  with  the  possible  exception  of  magnetohydrostatic  FBPS  which 
do  seem  to  be  stable. 

It  is  natural  to  try  and  explain  the  stability  characterisitcs  of  trial 
free  boundary  methods.  L.  Fox  [ 197  2]  suggested  the  following: 


Conjecture  2 


The  instabilities  of  the  trial  free  boundary  method  are  caused  by 
physical  instabilities  of  the  problems  being  solved. 


At  the  time  we  felt  that  this  conjecture  was  unsatisfactory  because  it 


relied  upon  the  rather  vague  concept  "physical  instability".  However,  the 
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evidence  of  instability  described  above  does  suggest  quite  strongly  that 
Conjecture  2 is  true. 

A fairly  decisive  test  for  Conjecture  2 is  provided  by  water-wave 
problems.  Southwell  and  Vaisey  | 1946,  p.  153]  describe  how  in  the  solution 
of  the  problem  of  the  flow  past  a planing  surface  the  computations  'took 
charge  and  a solitary  wave  emerged.  On  the  other  hand,  as  described 
above,  we  were  unable  to  solve  the  problem  of  a periodic  progressing  wave. 
This  suggests  that  solitary  waves  are  stable  and  periodic  progressing 
waves  are  unstable.  Recently,  Benjamin  and  Feir  [1965]  and  Benjamin  [1972] 
have  analysed  wave  waves  using  the  Korteweg-de  Vries  approximation  and 
have  come  to  the  same  conclusions  regarding  stability. 

An  important  consequence  of  Conjecture  2 is  that  since  viscosity 
usually  has  a stabilizing  effect  in  fluid  dynamics  we  are  led  to: 

Conjecture  3 

Viscous  fluid  FBPS  can  be  solved  using  trial  free  boundary  methods.  H 
The  importance  of  Conjecture  3 is  that  while  there  are  many  alternative 
methods  for  solving  inviscid  fluid  FBPS,  there  are  few  methods  for  solving 
viscous  fluid  FBPS.  The  recent  success  of  Tanner  [ 197  3]  and  Nickell, 

Tanner,  and  Caswell  [ 1974]  in  solving  viscous  jet  FBPS  using  the  trial  free 
boundary  method  offers  hope  that  Conjecture  3 may  be  true. 

We  conclude  with  some  remarks: 

(il  Southwell  and  Vaisey  [ 1946,  p.  122]  use  physical  arguments  to  explain 

the  instability  of  fluid  mechanics  FBPS  as  compared  to  porous  flow  FBPS. 
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(ii)  In  mechanics,  the  difference  between  stable  and  unstable  problems 
often  reveals  itself  in  the  governing  equations:  for  example,  the 


differential  equation  for  a damped  spring  contains  a damping  term. 
However,  for  many  FBPS  instability  seems  to  arise  through  in  a 
more  subtle  way  the  boundary  conditions:  the  governing  equations 
for  viscous,  stable  porous  flow  and  for  inviscid  unstable  fluid  flow 
are  both  Laplace's  equation.  One  can  of  course  surmise  that  the 
instability  arises  for  physical  reasons  such  as  the  introduction  of 
turbulence  at  the  boundary,  but  this  does  not  help  one  analyse 
the  problem. 
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4.  Concluding  remarks 

In  this  section  we  make  some  remarks  on  open  questions  about, 
and  possible  approaches  to,  trial  free  boundary  methods. 


4.1.  Error  estimates 

(k^  (k) 

Given  an  approximate  FB  I'  and  an  approximate  solution  u^ 

(k) 

one  would  like  to  estimate  the  error  in  the  solution,  u - u'  , and  also 

’ h ’ 

( k ) 

to  estimate  how  far  I is  from  F.  So  far  as  we  are  aware,  no  such 
error  estimates  are  known. 

There  is  a definite  need  for  error  estimates  because  even  experienced 

workers  find  it  difficult  to  estimate  the  error  in  the  solution.  A good 

example  of  this  is  provided  by  the  problem  of  an  inviscid  axially 

symmetric  wall  jet:  several  workers  (Trefftz  [ 1916]  , Southwell  and  Vaisey 

[1946]  , Rouse  and  Abul-Fetouh  [ 19  50]!  estimated  the  contraction 

coefficient  C to  be  .61  and  Rouse  [ 1962,  p.  188]  criticizes  Birkhoff 
c 

for  supporting  Garabedian's  estimate  = . 58  (Birkhoff  [1961,  p.  21]  1, 
but  the  best  current  estimate  is  = .59137  5 (Snider  [1971]). 

Some  authors  have  commented  on  the  fact  that  the  error  can  be 
quite  large  even  when  the  boundary  conditions  are  satisfied  to  within  a 
small  tolerance: 

(i)  Arms  and  Gates  [ 19  57,  p.  7]  observe  that  for  the  problem  of  a plane 

Riabouchinsky  cavity  the  initial  guess  F^  was  taken  to  be  the  exact 
analytical  FB.  Finite  differences  were  used  with  1145 
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gridpoinls  (see  Young,  Gates,  et  nl.  [ 195  5,  p.  12]  ).  The  velocity 

on  the  1'B  should  be  constant,  but  the  velocity  computed  from  u^ 

h 

vat icd  by  more  that  b%  , and  the  average  computed  velocity 
was  almost  2%  higher  than  the  true  value. 

(ii)  B.  W.  Hunt  [ 1967,  p.  18]  , who  used  integral  equations  to  solve  the 
problem  of  an  axially  symmetric  wall  jet,  compares  the  results  for 
two  curves  r^.  The  coordinates  of  the  curves  differed  by  less 
than  3%  at  every  point,  but  the  computed  velocities  behaved  quite 
differently:  for  one  curve  the  computed  velocity  differed  from  the 
asymptotic  velocity  by  less  than  1.  f>%  , while  for  the  other  curve 
the  computed  velocity  differed  from  the  asymptotic  velocity  by  as 
much  as  108%.- 

The  cases  quoted  above  are  perhaps  exceptional  because  the  flows  are 
badly  behaved  near  the  point  of  separation,  and  because  for  the  axially 
symmetric  jet  problem  the  fluid  velocity  on  the  FB  is  not  known  a-priori. 
Nevertheless,  these  cases  serve  as  a warning  that  relatively  high  errors 
can  occur. 

In  a few  instances  approximate  solutions  have  been  computed  for 
problems  whose  exact  solutions  are  known: 

(i)  Southwell  and  Vaisey  [ 1946,  Figure  11]  compare  the  exact  solution 
for  flow  through  a plane  Borda  mouthpiece  with  the  approximate 
solution  obtained  using  finite  differences  with  about  700  points. 

The  true  and  approximate  FBS  differ  by  about  5%. 
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(ii)  Analytical  and  numerical  solutions  are  available  for  porous  flow 
through  a rectangular  dam.  (Cryer[1976,  section  3. 1.1.1).  For 


a dam  of  height  24,  width  16,  and  downstream  depth  4,  the  height 
of  the  downstream  point  of  intersection  is  12.7132  to  four  decimals, 
and  the  best  value  obtained  using  a trial  free  boundary  method 
is  12.7  5 - an  error  of  . 03%. 

(iii)  j.  c.  Baker,  Beard,  and  Young  [ 1964]  use  the  self-consistent  method 
to  compute  the  earth's  magnetopause  for  three  simple  problems 
where  the  solution  is  known.  The  computed  results  are  indistinguish- 
able from  the  exact  results  to  the  accuracy  of  a graph. 

A substantial  part  of  the  machinery  necessary  to  estimate  the  error 
u - u^k)  is  already  available.  For,  to  obtain  error  estimates,  we  must 
be  able  to  estimate  two  quantities: 

(k)  (k)  (k) 

(i)  The  difference  u - u^  where  u^  is  the  approximate  solution 
of  the  problem 

tfu(k)  = 0,  in  f>(k), 

(1) 

£u(k)  = 0,  on  a&(k)  . 

(k)  (k) 

(ii)  The  difference  u - u where  u satisfies  (1)  and  u is 
the  solution  of  the  problem 

£7u  = 0,  in  & , 

(2) 

6 u = 0,  on  9&  . 


We  now  consider  these  questions  in  somewhat  greater  detail. 
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Estimation  of  u 


(k)  (k) 

" uh 
( k ) 

If  is  computed  oy  any  standard  method  such  as  the  method  of  finite 

h 

differences  or  the  method  of  finite  elements  then  there  is  an  extensive  literature 

( k ) ( k ) 

devofed  to  estimating  the  error  u ' - u,1  . There  are  often  several  different 

n 

approaches . 

For  finite  difference  methods  the  following  classical  approaches  were 


used:  Courant,  Friedrichs,  and  Lewy  [1928]  (variational  methods);  Gerscngonn 
[19  30]  (maximum  principle);  Petrovskii  [1941],  Jamet  and  Parter  [ 1967]  (barrier 
functions).  In  addition  there  have  been  many  other  approaches  to  finite  differ- 
ence methods:  Motzkin  and  Wasow  [19  53  ],  Bramble  and  Hubbard  [ 1964], 
Bramble,  Hubbard,  and  Thomee  [1969]  (monotone-type  approximations);  Nitsche 
and  Nitsche  [I960];  Cea  [1964],  Aubin  [ 1962]  (external  approximations); 

Lebaud  [ 1969],  Frehse  [1969]  (nonlinear  equations);  Lebaud  and  Raviart  [ 1969] 
(discontinuous  coefficients);  Schaeffer  [1972].  Diaz  and  Roberts  [ 19  52]  draw 
attention  to  the  similarity  between  iterative  methods  for  continuous  and 
discrete  problems. 

For  finite  element  methods  the  standard  approach  is  to  use  variational 
methods,  but  Ciarlet  and  Raviart  [to  appear]  have  used  the  maximum  principle. 

We  have  mentioned  all  the  above  methods  because  at  the  present  time 
it  is  not  clear  which  approach  will  be  best  for  FBPS. 

(k) 

Much  of  the  literature  is  based  on  the  assumption  that  d&  is 

smooth  and  is  therefore  not  always  applicable  to  FBPS  which  usually  involve 

(k) 

corners,  but  the  case  when  d&'  has  corners  has  been  considered  and 
references  are  given  in  section  4.  2. 

L 
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There  is  a considerable,  but  scattered,  literature  on  the  properties 


of  elliptic  difference  equations,  and  some  of  these  properties  probably 
play  a role  in  trial  free  boundary  methods.  See  Heilbronn  [ 1949]  , 

Kemeny  and  Snell  [ 1961]  , Varga  [ 1966]  , Ciarlet  [ 1970a]  . 

_ . , (k) 

Lstimation  ot  u - u 

(k) 

The  estimation  of  the  difference  u - u , that  is,  the  estimation  of  the 
variation  in  the  solution  of  a boundary  value  problem  with  respect  to  variations 
of  the  boundary  has  been  considered  from  several  viewpoints  in  the  literature: 

(i)  The  classical  variational  formulas  of  Hadamard  for  the  Green's 
function  (Garabedian  [ 1964,  p.  558]). 

(ii)  Estimates  for  the  variation  in  the  solution  of  the  Dirichlet  problem 

for  linear  elliptic  equations  of  order  m with  respect  to  variations 

of  the  boundary  have  been  obtained  by  Saak  [1972]  . 

(k ) 

(iii)  Since  £ and  £ differ,  it  is  desirable  (and  perhaps  essential) 

(k)  (k) 

to  be  able  to  continus  u and  u across  F and  F , respec- 
tively, so  as  to  be  able  to  compare  u and  u^1  on  the  same 
domain.  There  is  a growing  literature  on  numerical  analytical 
continuation:  G.  L.  Lewis  [ 1965],  Henrici  [1966],  K.  Miller  [ 1970]  . 

(iv)  The  question  of  domain  variations  arises  in  the  theory  of  finite 

elements,  because  in  general  the  boundaries  of  the  finite  elements 

do  not  coincide  with  the  boundary  of  the  domain  of  the  problem 

being  solved.  Berger,  Scott,  and  Strang  [ 1972]  , Strang  and  Berger 

[ 1974]  and  Thomee  [ 1973,  1973a]  give  estimates  for  the  differences 
( k ) ( \ 

u - u and  grad  (u  - u ) for  Poisson's  equation  in  the  plane. 

See  also  Necas  [1967,  p.  162). 
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(V) 


Steffen  [ 1972]  considers  FBPS  for  harmonic  functions  in  the  plane. 
Steffen  considers  the  boundary  conditions  on  r to  be  of  the  form 


Ftu 


du 

On 


jS 

*2  = 0. 


(vi) 


where  <>“  and  are  given  functions  which  depend  upon  x 
and  S,  . Steffen  shows  that  if  is  close  to  & then 

and  arc  almost  zero. 

Muhlig  [ 1073]  discusses  the  dependence  of  the  solution  upon  the 
domain  for  axially  symmetric  jets  from  curved  orifices. 


4.2.  Treatment  of  singularities 

Most  solutions  of  FBPS  have  singularities  because  most  FBPS 
involve  corners.  Here,  we  make  some  remarks  about  the  analytical  and 
numerical  treatment  of  singularities.  We  believe  that  this  is  an  area 
which  deserves  study  because  of  the  possibility  of  substantially  increasing 
the  accuracy  of  numerical  methods  for  FBPS. 

As  regards  the  analytical  study  of  singularities,  two  approaches 
have  been  used: 

(a)  Asymptotic  expansions 

Lewy  [ 19  50 , 19  57]  , Wasow  [ 19  57a]  , Laasonen  [ 19  57]  , Lehman  [ 19  59]  , 
Wigley  [1964,  1970,  1974]  have  derived  asymptotic  expansions  for  the 
solution  of  second  order  elliptic  equations  in  the  neighbourhood  of  corners. 
Kondrat’ev  [ 1967,  1970]  gives  expansions  for  general  equations  (of 


arbitrary  order) . 
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(b\  IntegrabiliU 


It  is  of  importance  to  know  in  which  Sobolev  space  the  solution 
lies.  See  Birman  and  Skvorcov  [ 1962]  , Kondrat'ev  [ 1967]  , Hanna  and 
Smith  [ 1967]  , and  Grisvard  [ 1971,  197  2]  , Kudrjavcev  [ 1974]  . 

See  also  Kellogg  [ 1966,  1970,  1971]. 

The  approximate  treatment  of  singularities  is  not  well  developed. 

We  list  below  several  techniques  which  have  been  used: 

(a)  Elimination  of  a singularity  by  conformal  mapping 

If  the  governing  equation  is  Laplace's  equation  in  the  complex 
z-plane  and  if  it  is  known  that  b ~ (z  - z^)  at  a corner  z^  = x^  + iy^, 

Of 

then  the  singularity  can  be  removed  by  the  conformal  mapping  w = (z  - z^j  . 
This  technique  was  applied  by  Mason  and  Farkas  [1971,  p.  19;  1972]  in 
conjunction  with  a trial  free  boundary  method. 

(bi  Subtracting  off  a singularity 

If  either  the  solution  or  the  FB  is  being  approximated  by  a combination 
of  functions  then  the  functions  should  obviously  be  chosen  so  as  to  match 
any  known  singularities.  Examples  will  be  found  in  Garabedian  [19  56; 
1956a],  Cryer  [ 1968 , p.  40],  Mason  and  Farkas  [1971,  p.  22],  J.  M.  Taylor 
[ 1971]  , Aitchison  [ 197  2]  , Fox  and  Sankar  [ 1969]  , Fox  [1971].  See  also 
Wigley  [ 1969]  . 

(c)  Construction  of  special  finite  difference  equations 

Motz  [ 1946]  , Thom  and  Apelt  [ 1961  j , L.  C.  Woods  [ 1953] . 
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w 


(d)  Dual  series 

Whiteman  [ 1967 ] , Whiteman  [ 1969 ] . 

(e  I Mesh  refinement  near  singularity 

Recently,  Volkov  [ 1966,  1966a,  1968]  has  analysed  the  error  for 
composite  meshes  consisting  of  a square  grid  of  gridlength  h "patched" 
to  'polar  grids  1 in  the  neighbourhood  of  corners,  but  we  are  not  aware 
of  any  computational  results. 

Proofs  of  convergence  for  finite  difference  methods  have  been  given 
by:  Walsh,  and  Young  [ 19  53 , 1954,  1957],  Wasow[1957],  Laasonen  [ 19  57a , 
1958],  Wigley  [ 1965]  , Kellogg  [ 1966]  , Bramble,  Hubbard,  and  Zlamal  [ 1968]  , 
Veidinger  [ 1968]  , Whiteman  and  Webb  [ 1970]  , Sultanova  [1971]  , Babuska 
and  Kellogg  [ 197  2]  . 

Finite  element  methods  for  problems  with  singularities  have  been 
considered  by:  Babuska  [1970,  1970a,  1974],  Babuska  and  Rosenzweig  [ 1972]  , 
Birkhoff  [ 197  2]  , Strang  and  Fix  [ 197  3,  chapter  8]  . 

Symm  [ 1973]  has  considered  integral  equation  methods  for  problems 
with  singularities. 

Almost  all  of  the  references  quoted  above  consider  fixed  boundary 
problems.  There  is  of  course  a great  need  to  consider  the  singularities 
at  the  points  of  detachment  of  the  FB.  This  will  be  discussed  in  detail 
in  a subsequent  survey.  Here  we  merely  mention  the  work  of  Aitchison  [ 197  2] 
(see  also  J.  M.  Taylor  [ 1971] ) who  obtains  an  expansion  of  the  FB  for  the 
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porous  fl  w through  an  earth  darn  in  the  neignbourhood  of  the  point  B 
at  which  the  FB  intersects  the  seepage  surface,  Aitchinson  uses  this 
result  in  two  ways:  (a)  to  help  move  the  boundary  (b)  to  compute  the 
solution  c near  B. 
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which  the  boundary  is  found  by  guessing  the  position  of  the  boundary,  solving  a 
boundary-value  problem  in  the  resulting  region,  and  using  this  solution  to  find  a 
new  guess  for  the  position  of  the  boundary,  the  procedure  being  repeated  until 

Tins  report  gives  a comprehensive 
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